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Section  I:  Introduction 

In  this  paper  Eric  Reissner's  equations  [8] 
for  axisymmetrical  deformations  of  thin  shells  of  revo- 
lution are  specialized  to  the  case  of  deformation  of 
a thin  ellipsoidal  shell  under  uniform  normal  pressure. 

i 

The  linked  pair  of  non-linear  differential  equations 
which  results  is  solved  approximately  by  the  Bubnov- 
Galerkin  method,  producing  the  first  complete  descrip- 
tion, including  collapse  and  deep  buckling,  of  an 
ellipsoidal  shell  under  pressure. 

An  oblate  ellipsoidal  shell  is  studied  because 
deformation  of  such  a shell  is  likely  to  occur  axisym- 
metrically.  Oblate  shells  under  pressure  have  aroused 
the  attention  of  Clark  and  Reissner  [lj  , who  deter- 
mined the  range  of  usefulness  for  such  shells  of  a 
linear  approximation  to  Reissner's  equations,  and  of 
Danielson  [2]  , who  used  special  buckling  equations 

to  determine  buckling  pressures  for  a wide  range  of 
geometries.  The  present  paper  extends  these  results 
through  a unified  and  refinable  discussion  of  all  states 
of  a single  shell.  The  comparable  numerical  results 
are  in  reasonable  agreement  with  Danielson's,  and  a new 
number  is  determined  — the  lower  critical  pressure. 
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the  minimum  pressure  which  can  sustain  the  shell  in  a 
buckled  shape. 

Three  technical  novelties  with  possible  applicabil- 
ity to  other  problems  are  worth  emphasis.  In  the  Bubnov- 
Galerkin  method  as  implemented  here,  sine  series  rather 
than  appropriate  eigenfunction  series  are  substituted 
for  the  dependent  variables:  in  effect  the  vanishing  of 

some  integrals  is  abandoned  in  favor  of  obtaining  all 
integrals  simply.  A powerful  addition  to  the  method  of 
continuation,  due  to  Rauch  in  [[6  3 and  useful  for 
taking  curves  around  corners  and  through  loops,  is  illus- 
trated here  as  well.  And  Polak's  stable  Newton-secant 
algorithm  for  solving  systems  of  equations  [43  is 
implemented  here  and  in  Rauch  et  al  [63  : tested 

FORTRAN  programs  used  in  both  works  are  included  in  the 
appendix. 


Section  II:  Establishment  of  Equations 

We  begin  with  the  equations  for  axisymmetric  defor- 
mations of  thin  shells  of  revolution  developed  by  Eric 
Reissner  in  [ 7]  and  ]_*8]  . 

In  these  equations*  (ro»^b,2(P  811,6  cylindrical 

coordinates  of  a generic  material  point  X on  the  central 
surface  of  an  undeformed  thin  shell  of  revolution;  ^ is 
the  angular  deviation  from  the  horizontal  of  the  radial 
tangent  to  the  central  surface  at  X;  $ is  a parameter 
upon  which  rQ  and  zQ  depend,  and  with  respect  to  which 
all  derivatives  are  taken;  ct  Q = V"^o  •>2  ♦ (v>2  is 
the  ratio  of  an  infinitesimal  length  of  meridian  to  an 
infinitesimal  change  in  S near  X. 

(r,^Q,z)  are  the  new  coordinates  of  X under  defor- 
mation (Recall  that  X is  a material  point,  capable  of 
motion.)  , and  ^ is  the  angular  deviation  from  the 
horizontal  made  by  the  radial  tangent  to  the  deformed 
central  surface  at  X.  The  variables  u,  w,  and  are 
defined  as  r - rQ,  z - zQ%  and  fQ  - <f>  , respectively. 

E is  Young's  modulus  for  the  material  of  the  shell, 
the  ratio  of  uniaxial  stress  to  the  extension  it  produces; 
if  is  Poisson's  ratio  for  this  material,  the  ratio  of 
transverse  contraction  to  the  extension  produced  by  uni- 
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axial  stress;  h is  the  thickness  of  the  shell;  and  C 
and  D are  defined  to  be  Eh  and  EhV(12  j]l-V2J  ). 

Py  and  are  vertical  and  horizontal  components 
of  the  pressure  on  the  shell. 

Finally,  consider  surface  tractions  across  infinites- 
imal surfaces  normal  to  the  radial  tangent  to  the  central 
surface;  let  these  be  integrated  with  respect  to  length 
along  a material  line  through  X and  perpendicular  to  the 
central  surface;  let  the  resultant  vector  be  resolved 
into  a vertical  component  V and  a horizontal  component  H 
(Rotational  symmetry  of  shell  and  deformation  are 
assumed,  and  so  no  other  component  is  present.);  then  ^ 
is  defined  to  be  *qH. 

With  these  definitions  in  mind,  we  discard  all  but 
the  three  most  important  terms  of  second  order  or  higher 
in  ^ and  ¥ from  the  equations  (III)  and  (IV)  of  ^8^  . 
We  also  discard  one  linear  term  involving  V,  and 

p 

rQ  p^.  We  assume  E,  h,  and  V,  and  hence  C and  D,  to 
be  constant  for  the  shell,  and  we  multiply  through  by 
rc/o<0*  This  leaves 


(1)  % + & )'f 


• . 


ro'  V 


Our  present  concern  is  with  an  ellipsoidal  shell 

of  revolution  whose  central  surface  is  given  in  cylin- 

2 2 2 2 

drical  coordinates  by  rQ  + zQ  /b  - a . When 


b * 1 the  ellipsoid  is  a sphere,  and  we  desire  equations 
(1)  and  (2)  to  take  the  form  of  Rauch's  (5a)  and 

(5b). 

For  simplicity  we  shall  choose  our  unit  of  length 


I 
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axial  stress;  h is  the  thickness  of  the  shell;  and  C 
and  D are  defined  to  be  Eh  and  Eh^/(1 2 £l-V^J  )• 

Py  and  pg  are  vertical  and  horizontal  components 
of  the  pressure  on  the  shell. 

Finally,  consider  surface  tractions  across  infinites- 

imal  surfaces  normal  to  the  radial  tangent  to  the  central 

l 

surface;  let  these  be  integrated  with  respect  to  length 
along  a material  line  through  X and  perpendicular  to  the 
central  surface;  let  the  resultant  vector  be  resolved 
into  a vertical  component  V and  a horizontal  component  H 
(Rotational  symmetry  of  shell  and  deformation  are 
assumed,  and  so  no  other  component  is  present.);  then  ^ 
is  defined  to  be  rnH. 

i 

With  these  definitions  in  mind,  we  discard  all  but 
the  three  most  important  terms  of  second  order  or  higher 
in  ft  and  from  the  equations  (III)  and  (IV)  of  [^8  J . 

We  also  discard  one  linear  term  involving  V , and 
r()2pH*  We  assume  E,  h,  and  V,  and  hence  C and  D,  to 
be  constant  for  the  shell,  and  we  multiply  through  by 
r(/o<0*  This  leaves 
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so  that  a ■ 1 and  then  choose  our  unit  of  force  so  that 


E - 1 

. Then, 

with  the  parameter  y defined  to  be 

(This 

parametrization  differs  from  that  used  in  £1]  • )» 

the  variables 

appearing  in  (1)  and  (2)  are  given  by 

(3) 

C * 

Eh  . h 

<*) 

T»  _ 

Eh3  h3 

12(1  -t/2)  12C1-V2) 

(5) 

tan  % 

dz0  _ _ 

*0  Z0 

(6) 

P_  a 

sin  £ 

o 

b2cos2  ? + sin2  5 

(7) 

2 ■ 

p 

-b  cos  £ 

zo 

Vb2cos2  ? + sin2f 

(8) 

V - 

p 

b C08$ 

(b2cos2*  + sin2?  )3^2 

(9) 

V - 

b2sin£ 

(b2cos2f  + sin2?  )3/2 

(10) 

*0  -V 

(b2cos2C  ♦ sin2?”  )3/f2 

(ID 

Is 
* 0 

-tt  sinC  (b2cos2f  ♦ sin2f  ) 
b2 

m 

7 

— — 


(12) 

■ 

cl  o 

2 

cos?  + 3 1 "*b-  sin2?  cos? 

b^ 

(13) 

o jo 

■ 

cos  £ 

(1<0 

'37* 

o jo 

B 

- sin  f 

(15) 

V 

*0 

sin  ? 

(16) 

o jo 

a 

o 

o 

CO 

(17) 

r ' 

ro 

b2  cot? 

ro 

b2cos2?  + sin2? 

(18) 

v 

b2 

ro 

b2cos2?  + sin2? 

(19) 

%\  4 

o lo 

°l°. 

b2cos2? 

sin  ? (b2cos2  f ♦ sin2  f ) 

(20) 

ro  ro'zo 

i p * 

b coss 

wo  TT 

® 5 5 ^ 

b^cos^  + sin*1? 

and 
PH  - 

(21) 


Under  uniform  normal  pressure  />  , p ^ sin 

Py  - />  cos  • We  use  the  approximations 
- f sin  ? and  py  - p cos  f and  obtain 


r02*H 


b2cos2  £ + sin2  ? 


- — 


(22) 


^r02pH^' 

-/?(3b2sin2/  cos2?  + (1  ♦ 2b2)sin^t  cosf 
(b2cos2/  + sin2/  )2 


(23) 

<r0V)- 

(24) 

(r0T) 

-r^« 


o'cto 


p >. 

-/^b  sin  T cos  £ 
(b2cos^  T + sin  ( )‘ 


2(1  - b2) (b2cos2  / + sin2 i ) 2(1 

-/’sin2/ 

■ ■ ■ — 3 — tp—  r 3 

2 ( b2cos2  / + sin2  / ) 


The  peculiar  integration  constant  -1/(2(1  - b2))  was 
introduced  so  that  all  the  variables  in  (3)  - (24) 
approach  the  corresponding  variables  in  Rauch  £5  J 
as  b approaches  1. 

Substitution  of  (3)  - (24)  into  (1)  and  (2)  produces 


(25) 


r /S" 


sin2f  cos  / p ' 


^sin^  + — sin^rj^" 

+ j^cos£  + 5 1 ~^b  sin2f  cos  ^ ' 

. f ll + J .ini]  / 

j^sin  / (b2cos2  ? ♦ sin2/)  J > 

» 12(1  - V2)  _ b2sint  V 

b5  (b^cos^  / ♦ sin2  X ) 


b2 

.2 2 <■ 


gb  cf  ? j—  ♦ V sin{ 

sin  f (b2cos2?  ♦ sin2/  ) 

12(1  - V2)  f b2sinf  ^ 


V 


f. 


(26) 


fi  b2sin2f  cos? 
2(b2cos2f  + sin2J  'y>> 

/>  b2sin5f  & 

'5  ■■  «■  5 ■ w' 

2(b^cos^f  + sin^f  )• 
i - b2 


b2cosf  @ ¥ 
(b2cos2£  + sin2C  )‘ 


Vjy  • ' 


cos  f + 3 


1 - bc  .2 


sin^T  cos  f vj/' 


b2cos2  f 

sin  X (b2cos2  £ + sin2,  ) 

h I b2sinj>  $ 

(b2cos^J  + sin2  f )• 


•— 1 

5J  f 


2(b<icos<:? 


b2cos?  6 2 1 

‘ a(b^»  ♦ 

<7wTf~S7F7  [-  * 

♦ b2(3  - 2i/)8in2|’  cos^T 
+ (1  - j V + 2b2)sin*T  cos,] 


If  we  let  b ■ 1 in  equations  (23)  and  (26), 
divide  through  by  sin  X , and  sake  the  substitution 
/'a  2^1  cos  f sin  X ♦ ^ * we  obtain  the  equations 
(3a)  and  (3b)  of  £3]  as  desired. 


». 


allow 


Reissner's  formulas  (10)  and  (12)  of  [8^ 
the  determination  of  the  horizontal  and  vertical 
displacements  u and  wt  once  fl%  ^ % 4/*  ¥'»  and  (rQV) 
are  known.  In  the  present  case  the  formulas  are 

(27)  u - j;  ^ - *q2/°  ain  / - V^cos<^ 

- V (r0V)sin^ 
f r «<nsin  9^  C 

(28)  w - J 1^— £r  h~  | *0h  * Xcos  <f>  + (r0V)sin 

- V^'  - rQ2V  f sin  - *Qsin  df 

The  integral  for  w is  normalized  so  that  the  south  pole 
shows  no  vertical  displacement;  this  integral  can  be 
approximated  adequately  by  the  trapezoidal  rule. 


Section  III:  Algebraization  of  the  Equations 


We  proceed  to  find  approximations  toft  and  4^ 

through  the  Bubnov-Galerkin  method:  we  replace  all 

the  dependent  variables  and  their  derivatives  by 

m 


truncated  series,  ft 


B j 8 in  d X ♦ 


' - d B.  cos  dt  » 


X.  (-d^)  B.  sin  d?  » 


^ ■ 51  Pj  sin  d F < 

J d*i  0 


••  • 


d Pd  cos  d 5 , 


51  (“d2)  Pi  sin  d F . 

d-i  3 


This  leaves  ub  two  equations  involving  only  the 
parameters  ht/f,  b and  i/ , the  independent  variable  | , 
and  the  2m  unknown  constants  B^  and  P^«  Ve  eliminate 
T by  multiplying  each  equation  by  sin  if  and  integrat- 
ing both  sides  from  0 to  ft • If  this  is  done  for  i 
from  1 to  m,  we  are  left  with  2m  algebraic  equations 
involving  only  the  2m  unknowns  Bj  and  Pj , the 


parameters,  and  a number  of  evaluated  integrals.  It 
is  this  system  of  2m  equations,  rather  than  (25)  and 
(26),  which  we  attempt  to  solve. 

Ve  write  this  system  as  equations  (44^)  and  (45^) 
in  terms  of  the  integrals  and  their  combinations 

which  we  define  below. 


(29)  S1(i, j) 


(30)  S2(i,d) 


(3D  S3(i,d) 


(32)  S^(i,d) 


(33)  S5(i,j) 


sin?  sin  if  sin  j?  d? 


Isin^?  sin  i£  sin  d5  d? 
0 


cos?  sin  i?  cos  d?  d? 

'0 

C?  2 

j sin  ^ cos?  sin  i?  cos  dT  d'S 
10 

cos  T 

q sin?  (b2cos2?  ♦ sin2?  ) 

sin  iS  sin  d5  d? 


(34)  S6(i,d)  - f 7-5 2 sin^ — 2 — r: 

b Jn  (b2cos2S  + sin2?  )' 


sin  if  sin  d$  d? 


(35)  S?(i) 


sin  <,  cos 


(b^cos^J  * sin2?  )' 


sin  iX  &t 


(36)  S8(i,d,k) 


(b2cos2T  ♦ sin2S  )' 


(37)  S9(i,j) 


sin  iS  sin  sin  k£  d< 

' 3 
sin^ 

(b2cos2*s  + sin2$ 


(38) 


(39) 


S10(i) 


sin  i5  sin  &% 

p 

sin  % cos$ 
(b2cos2iS  + sin2$  )‘ 


f'*'  ? 

Su(i)  . rSifi.  Si 

11  Jr\  (b2cos2?  •*> 


siu  5 Y 


(40)  s12(i) 


sin^T  cos? 
(b2cos2?  + sin2?  )‘ 


(41)  01(i,d) 


1 - b‘ 


(42)  C2(i,J)  - C1(i,J)  ♦ 2Vs,(ifJ) 


(43)  C5(i)  - - £ b2  S10(i)  + b2  (3  - £V)  Sn(i) 

+ (1  - ^ * 2b2)  S12(i) 

Equations  (25)  and  (26)  now  give  rise  to 

(4V  - ' 2il ^2)b2  [-^s7(i) 

+ ZL  £ - S6(i,d)  Pj  - £jk/>S9(i,d) 

* s8(i»j»k)  pk]  1 

(45;:)  Z1  C2(i,o)  P,  - hb2  I 5”  js6(i,d) 

d-i  ^ J Ld-i  t b 

*2  £ s8(i-l-k>  Sj]  + /°3(i) 

Once  the  and  are  known  for  a particular  pair 
of  values  of  b and  v,  equations  (45^)  express  the  P^ 
in  terms  of  the  B^.  Thus*  solving  this  system  reduces 
to  finding  the  m variables  B^.  We  proceed  to  develop 
an  inexpensive  way  to  compute  the  and  C^. 
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Section  IV:  Evaluation  of  Integrals 


Except  for  Sr,  the  integrals  S,  are  all  of  the  form 
f'i'  ? w 

(Kernel) (Product  of  an  odd  number  of  sines)  /\ 

J°  \y 

(Product  of  0 or  more  cosines)  d % . 

Any  of  these  may  be  evaluated  with  reasonable  economy 

i 

by  a two-step  process:  (1)  use  Simpson's  Rule  to 

f •* 

evaluate  (Kernel)  sin  1$  dl  for  a sufficiently 

Jo 

wide  range  of  values  of  the  integer  I;  (2)  express  the 


S^  in  question  as  a simple  function  of  the  integrals 
so  evaluated  through  use  of  the  formula 

(46) 


sin  1^  sin  I2  . • . sin  In  cos  cos  J2  • • • cos 

/ •»  \ £n/2^  r"" ’ 

• ('2nU-l  • - fWV-  • -ljt> 

In  this  formula  (n/2]|  is  tho  greatest  integer  in  n/2; 
"trig”  is  "sine"  if  n is  odd,  "cosine”  if  n is  even; 
the  sum  is  taken  over  all  possible  combinations 

of  plusses  and  minusses  in  the  argument  of  trig;  and 
the  sign  of  trig  is  the  product  of  the  signs  of  1^,  I2, 


. . • In  in  its  argument. 


For  example,  to  evaluate  Sg(i,j,k),  first  evaluate 
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sin  I?  d.C  /(b2cos2t  + sin2?  )^2  for 


T(f)  - j 

all  integers  I from  2 - 2m  to  3m  + 1.  Then  Sg(i,j,k)  « 
- g £ T(i+j+k+l)  + T(i+j+k-l)  - T(i+j-k+l)  - T(i+j-k-l) 

- T(i-j+k+l)  - T(i-j+k-l)  + T(i-;j-k+l)  + T(i-j-k-l)^ 
The  economy  of  this  process  is  increased  when  the 
symmetries  and  anti-symmetries  of  the  kernel  functions 
and  of  sin  1$  are  considered.  All  of  the  kernels  and 
all  odd  sines  are  symmetric  about  1?/2;  so  for  I odd, 
Simpson's  Rule  need  only  be  applied  to  half  the  interval 
0 toTl^  . On  the  other  hand,  even  sines  are  anti- 
symmetric about  <lt'/2,  so  that  for  I even, 

(Kernel)  sin  1$  d|  * 0,  and  Simpson's  Rule  need 

not  be  invoked*  Finally,  sin  I?  « - sin  (-1)?  so  that 

Simpson's  Rule  need  only  be  used  on  integrals  of  the 
r1?/l 

form  ] (Kernel)  sin  IX  d?  where  I is  a positive 

J 0 

odd  integer. 

^ The  process  described  above  fails  for  S^(i,j)  * 

cos2?  sin  it  sin  d?  /(sinT  (b2cos2$  + sin2*  )), 

,J0 

since  for  this  integral,  formula  (46)  would  suggest 
an  evaluation  of  integrals  of  the  form 


I 


»r 


'o 

step 


cos  ltd?  /(sin?  (b2cos2?  + sin2?  ))  as  a first 


Unfortunately,  some  of  these  integrals  do  not 
converge,  e.g.  b ■ 1,  I ■ 1.  So  the  integrals 
are  evaluated  by  Simpson's  Rule  directly.  Time  is 
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saved  however  when  symmetry  considerations  lead  to  the 
observations:  (1)  S5(i,J)  . S5U,i);  and  (2)  when 

i + j is  odd,  S^(i,j)  ■ 0.  In  computing  the  integrals 
S^,  use  is  made  of  the  formula 

sin  IT  sin  (I  - 1) C 

: " ■ X * cos  J + COS  (1-1)5 

sin  5 sin  S 
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Section  V:  Polak's  Method 

In  this  section,  x is  a vector  whose  components  are 
(x-^Xg,  ...  xffl),  y is  a vector  whose  components  are 
(y1»y2»  ...  ym)»  etc. 

Newton's  method  for  solving  a system  of  m equations 
in  the  unknowns  x^,  Xj,  . • . xm  requires  that  each 
equation  be  put  in  the  form  fi(5c)  ■ 0.  If  the  guess  x 
is  near  an  unknown  solution  y,  Taylor's  Theorem, 
truncated  to  exclude  non-linear  terms,  gives 

(y1-x1)  + . . . + I (ym-*m)  o 

m|- 

Let  this  linear  system  be  solved  for  the  correction 
z » y - x,  and  let  z be  added  to  the  original  guess  x 
resulting  in  a new  vector  x*.  In  many  cases  x*  will 
be  significantly  closer  to  the  true  solution  y than  was 
x.  If  x*  is  used  as  a new  guess, and  if  this  process 
is  repeatedly  applied,  and  if  the  original  guess  was 
sufficiently  near  the  true  solution,  convergence  to  the 
true  solution  might  possibly  occur.  When  convergence 
does  occur,  it  is  usually  rapid. 

Polak  £4}  uses  a secant  approximation  to  the 
first  derivatives  required  by  Newton's  method,  and 
he  keeps  Newton's  method  from  violent  instability 
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by  requiring  that  for  any  now  guess  the  residual 

+ f2(x*)2  + . . . + fm(x*)2  be  smaller  than 
the  corresponding  residual  for  the  last  guess,  x — 
else  the  new  guess  is  not  accepted.  To  keep  the  method 
from  coming  to  a halt  in  such  a case,  Polak  takes 
advantage  of  residuals  that  might  as  well  be  computed 
in  the  course  of  computing  the  first  derivatives 
necessary  for  Newton's  method:  for  all  the  f^'s  must 

be  computed  at  (x1,x2,  ...  x^,  ...  xm) 

+ (0,  0,  ...€,...0)  for  each  k, 

and  it  is  possible  that  such  a vector  might  be  an 
improvement  over  the  guess  x. 

Further  details  and  a FORTRAN  program  for  the 
implementation  of  Polak 's  method  are  in  the  Appendix. 
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Section  VI:  Continuation 


r 


Our  system  of  equations  involves  the  parameters 
b,y0,  h,  V , and.  mt  and  the  variables  B^.  Suppose 
that  for  a particular  assignment  of  values  to  the 
parameters,  a solution  for  the  is  known.  Then  let 
us  take  this  solution  as  our  initial  guess  for  a new 
problem  in  which  the  assignment  of  values  to  one  of 
the  parameters,  say^o  , is  only  slightly  different  from 
the  assignment  in  the  last  problem,  and  otherwise  the 
assignments  are  the  same.  If  we  do  so,  and  apply  Polak's 
method,  our  chance  of  finding  a solution  to  the  new 
problem  is  good.  We  say  we  are  performing  a 
continuation  of  the  first  type,  with  p varying. 

There  is  a new  and  useful  continuation  of  a second 
type,  in  which  one  of  the  variables,  say  B^,  switches 
roles  with  one  of  the  parameters,  say  fi  • Suppose 
again  that  for  a particular  assignment  of  values  to 
the  parameters,  a solution  for  the  B^  is  known.  Then 
let  a new  problem  be  defined  by  assigning  the  same 
values  to  b,  h,  V , and  m,  varying  the  former  solution 
for  B^  slightly,  and  then  demanding  that  a set  of  values 
for  p , Bg,  Bj,  • • . Bm  be  found  to  satisfy  the 
equations. 
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The  second  type  of  continuation  has  been  extremely 
useful  in  several  situations.  The  graphs  in  Figures 
2,  3»  have  long,  nearly  horizontal  stretches  which 
would  have  been  exorbitantly  expensive  to  obtain  by 
doing  a continuation  of  the  first  type  with  p varying. 
Moreover,  the  loops  and  near-cusps  in  these  graphs 
would  have  been  impossible  to  discover  without  second- 
type  continuation.  Also,  transitions  from  even  (solid- 
line) solutions,  in  which  all  the  odd-numbered  are 
zero,  to  odd  (dotted-line)  solutions,  in  which  some  of 
the  odd-numbered  are  not  zero,  were  effected  by 
continuations  of  the  second  type,  with  an  odd-numbered 
varying.  (Odd  and  even  solutions  are  discussed 
further  in  Section  VII.)  And,  in  a related  problem, 
when  it  was  found  necessary  to  increase  m from  a point 
near  the  bottom  of  the  curve,  second-type  continuation 
with  a small  unimportant  B^  varying  hardly  at  all,  had 
to  be  combined  with  first-type  continuation  with  m 
varying,  because  of  the  sensitivity  of  the  system  to 
the  necessarily  integral  jumps  in  m. 

First  type  continuation  was  useful  in  an  interesting 
way  in  obtaining  one  of  the  two  starting  points  used  for 
the  development  of  the  curves  in  Figures  2,  3*  4. 

In  £5]|  and  £6}  attention  is  focussed  on  deformations 
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of  thin  spherical  shells,  and  variables  called  play 
roles  analogous  to  the  roles  of  the  in  the  present 
problem.  We  shall  call  the  modes  in  the  present 
discussion.  It  was  expected  that  in  some  solutions 
to  equations  (A1)  and  (B1)  of  £5}  » the  A^  would 

gradually  increase  as  i increased,  reaching  a maximum 
of  Aj  for  some  j,  then  gradually  decrease  to  a minimum 
of  A^  for  some  k,  then  gradually  increase  again  to  a 
relative  maximum  at  An  for  some  n,  etc.,  with  |A^| 
considerably  larger  than  |Ak|  which  itself  would  be 
considerably  larger  than  jAn|  , etc.  We  called  A^ 
the  critical  mode  in  the  solution,  and  remarked  that 
as  the  ratio  of  radius  to  thickness  ^ increased,  j would 
also  increase,  and  more  modes  would  be  necessary  to 
establish  an  acceptable  solution. 

It  was  found,  further,  that  (A')  and  (B*)  of  [5] 
could  be  solved  explicitly  when  the  number  of  modes 
was  2.  However  for  large  attempts  to  continue, 
varying  the  number  of  modes  from  2 up  to  and  past  j, 
met  with  consistent  failure.  It  seemed  that  the 
explicit  2-mode  solutions  would  be  useless  for  producing 
an  acceptable  solution  to  a reasonable  problem, 
involving  large  ^ — that  is,  involving  thin  shells. 
Before  abandoning  the  idea,  however,  we  decided  to 
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try  starting  with  an  unphysically  small  f ” ^ * for 
which  the  critical  mode  was  A^,  and  then  alternating 
first-type  continuations,  varying  ^ for  a while, 
then  varying  the  number  of  modes,  then  varying  ^ 
again,  never  allowing  ^ to  get  so  large  that  it 
required  a critical  mode  whose  subscript  was  larger 
than  the  current  number  of  modes.  This  worked,  and 
it  produced  the  solution  in  the  first  column  of  Table  1, 
which  is  listed  in  a form  translated  into  the  terms  of 
the  present  paper  in  the  second  column.  A refinement 
is  listed  in  the  third  column. 

The  solution  in  the  third  column  of  Table  1 was 
then  first-type  continued,  varying  m,  then  varying  b, 
then  varying  />  , then  varying  b again,  to  produce  a 
20-mode  starting-point  solution  for  an  ellipsoidal 
shell  1 unit  in  radius  at  the  equator,  1 unit  tall 
from  pole  to  pole,  .02  units  thick,  and  at  a 
pressure  of  .9E-4  units.  This  starting-point  was 
continued  to  produce  the  dotted-line  curves  in  Figures 
2,  3, 

The  solid-line  curves  were  all  continued  from  the 
trivial  solution  ■ 0 for  all  i to  the  problems 
p m o,  b • .5  , V ■ .3*  ■ ■ 20,  30,  or  40. 

The  solid-line  and  dotted-line  curves  all  have 
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common  solutions  at  their  highest  points.  These 
solutions  are  listed  in  Table  2,  and  may  be  continued 
to  produce  both  the  solid-line  and  the  dotted-line 
graphs . 

Numbers  are  listed  in  the  tables  in  the  computer 
version  of  scientific  notation:  in  this  notation 

•2910E-3  stands  for  .2910x10”^. 

It  should  be  remarked  that  the  close  agreement 
between  columns  2 and  3 of  Table  1 is  a good  verification 
of  the  consistency  of  the  present  paper  with  []6j  and  of 
the  correctness  of  the  computer  programs  used  in  both 
projects. 


1 

Table  1. 

Starting-point  solution 

Solution  to  (A') 
and  (B* ) of  [f\ 

Translation  to 
terms  of  the 
present  paper 
b - 1 

Refinement 
b - 1 

a 

E " 

50 

h - .02 

h - 

.02 

r 

■ 

.6 

f>  - .2910E-3 

/°  • 

■ .2910E-3 

Ai 

■ 

.2004E-3 

Bx  . . 1888E-2 

B1 

- .1879E-2 

A2 

a 

.3517E-3 

B2  - .3693E-2 

B2 

- .3676E-2 

A3 

a 

.5194E-3 

B?  . .5731E-2 

B3 

- .3707E-2 

a4 

a 

.7283E-3 

B4  - .7496E-2 

- .7463E-2 

A5 

m 

.97^7E-3 

B5  - .9716E-2 

B5 

- .9676E-2 

6 

a 

•1280E-2 

B6  - .1138E-1 

B6 

- .1133E-1 

1 A7 

a 

• 1628E-2 

B?  - .1362E-1 

B7 

- .1357E-1 

A8 

a 

.2015E-2 

B8  • .1480E-1 

B8 

- .1475E-1 

A9 

a 

.2377E-2 

B9  - .1651E-1 

B9 

- .1645E-1 

A10 

a 

. 2662E-2 

B10  - .1632E-1 

B10 

- .1623E-1 

A11 

a 

.2765E-2 

BX1  - .1655E-1 

B11 

- .1648E-1 

1 A12 

a 

.2672E-2 

B12  - .1410E-1 

B12 

« .1404E-1 

A13 

a 

.2393E-2 

B13  - .1279E-1 

B» 

- .1274E-1 

A14 

a 

.2022E-2 

B14  - .8461E-2 

BW 

- .8424E-2 

a 

.1624E-2 

B15  - .7036E-2 

- .7008E-2 

Residual  - 

j .7071E-4 

Residual  • 

.1407B-1 

Residual  - 

.5513E-4 
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Table  2.  Highest 

-point  solutions,  b » .5 

, h ■ 

• .02,  ■ .3 

m ■ 20 

m - 

30 

m * 

40 

f>  - .1443E-5 

/*  " * 

1572E-3 

/o  - .1524E-3 

b2  - .4389E-1 

B2 

m 

.4552E-1 

B2 

- .4590E-1 

- .2076E-1 

B4 

m 

.3203E-1 

B4 

- .2689E-1 

B6  - .1343E-1 

B6 

m 

.1978E-1 

B6 

- .1599E-1 

Bq  . .1003E-1 

B8 

m 

.5293E-2 

B8 

- .5709E-2 

B10  " *2973E-2 

B10 

m 

-.4477E-2 

B10 

- -. 1614E-2 

b12  - -.2934E-2 

B12 

m 

-.6866E-2 

B12 

- -.4335E-2 

B14  - -.5904E-2 

B14 

m 

-.3890E-2 

B14 

- -.3279E-2 

B16  - -.6002E-2 

B16 

m 

.9848E-3 

B16 

- -.3955E-3 

B18  " -*4287E“2 

B18 

m 

.5121E-2 

B18 

- .2587E-2 

B2q  - -.1957E-2 

B20 

u 

.7345E-2 

B20 

- .4699E-2 

B22 

m 

.7635E-2 

B22 

- .5678E-2 

B24 

m 

.6551E-2 

B24 

- .5692E-2 

B26 

m 

.4777E-2 

B26 

- .5075E-2 

B28 

m 

.2870E-2 

B28 

- .4154E-2 

B30 

m 

.1195E-2 

B30 

- .3172E-2 

B32 

- .2275E-2 

BM 

- .1528E-2 

B36 

- .9431E-3 

B38 

- .5105E-3 

B40 

- .1998E-3 

Odd-numbered 

are  0 and 

are  not  shown. 
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Figure  2.  Pressure  p versus  deflection 


OE  respectively 


Figure  3.  Pressure  p versus  deflection 


Pressure  fi  versus  deflection 


and  OE  respectively 


Section  VII:  Interpretation 


An  even  solution  to  this  system  of  equations  is 

m 

one  in  which  ft  • £ B.  sin  i ? is  anti-symmetric 

1?  i-1  1 

about  ^ . This  state  is  characterized  by  a symmetric 

alignment  of  the  shell  with  respect  to  the  equatorial 
plane,  and  by  the  vanishing  of  the  odd-numbered  B^. 

A solution  which  is  not  even  is  called  odd.  In  Figures 
2,  J,  4-  the  odd  solutions  are  indicated  by  dotted  lines, 
the  even  solutions  by  solid  ones. 

For  physical  reasons,  any  odd  solution  must  have  a 
matching  odd  solution  in  which  the  deformed  shell  looks 
like  the  original  deformed  shell  turned  upside-down. 
This  would  be  effected  by  a matching  solution  in  which 
the  even-numbered  are  the  same  as  in  the  original 
solution,  but  the  odd-numbered  B^  are  minus  what  they 
were  in  the  original  solution.  Any  odd  solution  and 
its  mate  produce  the  same  deflection  and  so  are  repre- 
sented by  the  same  point  in  Figure  2,  3»  or  4. 

The  graph  of  the  odd  solutions  may  be  expected 
to  have  an  even  solution  as  one  endpoint.  It  is  of 
great  interest  that  this  even  solution  happens  to  be 
the  highest  point  of  the  figure:  strings  of  odd 

solutions  have  been  continued  up  to  such  points  by 


I 
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continuations  of  the  second  type,  and  have  been  found 
to  transform  themselves  smoothly  either  into  strings 
of  even  solutions,  or  through  a single  even  solution 
into  strings  of  matching  odd  solutions  — depending 
upon  whether  the  particular  varying  was  even-  or 
odd-numbered.  Conversely,  a string  of  even  solutions 
for  the  problem  m - 20,  b » .6,  V > .3,  h - .02  has 
been  transformed  at  its  highest  point  into  a string  of 
odd  solutions  by  varying  an  odd-numbered  B^. 

Deflection,  shown  on  the  horizontal  axis  in 
Figures  2,  3,  4,  is  defined  in  this  paper  as  the 
difference  between  the  lengths  of  the  polar  axes  of 
the  undeformed  and  the  deformed  shells.  Since  the  poles 
cannot  pass  through  each  other,  deflection  cannot  be 
greater  than  the  length  of  the  undeformed  polar  axis, 
and  so  our  graphs  stop  at  this  length.  But  mathematical 
solutions  with  this  unphysical  characteristic  do  exist: 
the  graphs  seem  to  continue  to  rise  to  the  right  beyond 
the  realistic  cut-off. 

We  have  not  attempted  an  energy  analysis  of  our 
solutions.  But  it  seems  plausible  that  under  a given 
pressure,  a shell  with  a substantially  larger  deflection 
possesses  more  elastic  potential  energy  than  a shell  with 
a smaller  deflection.  Under  this  assumption,  the  shell 
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would  be  expected  to  follow  the  leftmost  accessible  branch 
of  such  a curve  as  the  one  in  Figure  4 when  in  an  unstable 
state.  With  this  tentative  principle  in  mind  we  offer 
an  interpretation  of  Figure  4. 

We  begin  at  the  origin  of  Figure  4,  corresponding 
to  the  first  meridian  of  Figure  5.  As  the  pressure  is 
raised  from  0 to  129,  the  deflection  increases  almost 
linearly  to  about  5/100  of  an  equatorial  radius.  If 
the  pressure  passes  this  crisis  point,  which  corresponds 
to  the  near-cusp  in  Figure  4 and  to  the  second  meridian 
of  Figure  5,  the  leftmost  available  solution  will  be 
on  the  left  side  of  the  breaking-wave-like  piece  of 
graph  near  the  top  of  the  figure.  If  now  the  pressure 
continues  to  increase  past  152,  which  corresponds  to 
the  highest  point  in  Figure  4 and  to  the  third  meridian 
of  Figure  5,  the  leftmost  available  solution  will  be  on 
the  dotted-line  graph  far  to  the  right  of  what  has  been 
drawn,  and  the  shell  will  collapse. 

On  the  other  hand,  if  the  pressure  passes  the  crisis 
point  and  then  is  reduced  slowly,  the  shell  might  find 
itself  on  the  near-vertical  solid-line  below  the  breaking 
wave.  By  the  time  the  pressure  had  decreased  past  80, 
the  shell's  poles  would  be  approaching  each  other  visibly, 
and  would  approach  each  other  very  dramatically  as  the 
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pressure  dropped  from  40  to  25:  the  deflection  would 

increase  from  14/100  to  48/100  in  that  range.  If  the 
pressure  then  dropped  below  25 » the  only  available 
solution  would  be  on  the  part  of  the  graph  near  the 
origin  — in  effect,  the  shell  would  snap  out  from 
the  shape  indicated  by  the  fourth  meridian  in  Figure  5 
to  a shape  close  to  the  undeformed  shape  of  the  first 
meridian. 

If  instability  is  a feature  of  the  solutions  on 
the  section  of  solid-line  graph  we  have  just  been 
discussing,  then  at  any  disturbance  the  shell  might 
leave  that  solid-line  graph,  with  its  symmetric  north- 
south  indentations,  and  seek  the  lower-potential 
further-left  dotted-line  graph  with  its  one-pole 
indentation.  Suppose  that  so:  then  as  the  pressure 

is  reduced  from  40  to  25,  the  deflection  would  increase 
from  8/100  to  26/100  before  snapping  out  from  the  shape 
shown  in  the  fifth  meridian  of  Figure  5* 

It  is  remarkable  that  the  bottom  halves  of  the 
fourth  and  fifth  meridians  of  Figure  5 match  so  well, 
and  that  in  Figure  4 the  even  graph  from  (14,40) 
rightwards  is  a half-speed  imitation  of  the  odd  graph 
from  (8,40)  rightwards.  (A  similar  comment  holds  for 
Figures  2 and  3 as  well,  and  for  graphs  obtained  in 
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Rauch  et  al  [6]  for  spherical  shells.)  It  is  remark- 
able also  that  the  snap-out  pressures  on  the  odd  and 
even  graphs  in  Figure  4 agree  to  three  places:  it  seems 
reasonable  to  take  either  number,  as  listed  in  Table  3, 
to  be  our  approximation  to  the  lower  critical  pressure, 
the  least  pressure  which  can  maintain  a buckled  shape 
of  the  shell. 

Danielson's  buckling  pressure  for  the  oblate 
ellipsoidal  shell  of  our  work  is  .137E-3  as  indicated 
in  Figure  10  of  ([2]  • This  corresponds  well  to  our 
crisis  pressure,  .1297E-3. 

Our  collapse  pressure  and  our  lower  critical 
pressure  are  118#  and  19$  respectively  of  our  crisis 
pressure.  The  first  onset  of  substantial  non-convexity 
occurs  as  the  pressure  decreases  (after  having  passed 
the  crisis  pressure)  in  the  range  71#  - 67#  of  the 
crisis  pressure. 

The  programs  which  performed  the  work  of  this 
paper  were  run  on  an  IBM  370/168  at  the  City  University 
of  New  York* 8 University  Computer  Center.  The  even 
m ■ 40  graph  was  obtained  using  100K  and  30  minutes; 
the  odd  m » 40  graph  required  340K  and  170  minutes. 

If  further  work  shows  that  the  most  important  pressures 
— the  crisis  pressure,  the  collapse  pressure,  and 
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Table  3«  Important  Pressures: 


b - 

.5,  h - .02, 

V - .3 

m a 20 

m ■ 30 

m ■ 40 

Crisis 

.1302E-3 

.1300E-3 

.1297E-3 

Collapse 

.1443E-3 

.1572E-3 

.1524E-3 

Outward-snap 

(even) 

.2576E-4 

.2485E-4 

•2488E-4 

Outward-snap  (odd) 

.2470E-4 

.2478E-4 

.2481E-4 

These  pressures  are  given  in  the  units  of  the  paper. 

To  convert  to  real  units,  multiply  by  Young's  modulus 
for  the  material  in  question.  For  example,  if  a shell 
fitting  the  description  b « .5,  h - .02  is  made  of 
a material  for  which  V » .3  end  E » 30,000,000  pounds 
per  square  inch,  then  such  a shell  would  collapse  at 
a pressure  near  4572  pounds  per  square  inch,  and  would 
be  incapable  of  sustaining  a buckled  shape  at  pressures 
below  744  pounds  per  square  inch. 
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snap-out  pressure  — all  continue  to  be  obtainable 
from  the  even  graphs  alone,  then  future  work  on  the 
question  of  deformation  of  ellipsoidal  shells  may 
concentrate  on  developing  only  these  much  more  economical 
graphs.  Future  work  should  also  attempt  to  include 
more  of  the  terms  of  £•  Reissner's  differential 
equations. 


Appendix:  FORTRAN  Programs 


I 

I 

This  appendix  serves  primarily  as  a place  of  deposit 
for  the  FORTRAN  programs  which  developed  the  data  used 
in  the  present  paper  and  in  Rauch  et  al  [6j  • These 
programs  fall  into  two  classes:  major  programs  which 
utilize  a particular  subprogram  SECANT  to  solve  systems 

I of  equations;  and  satellite  programs  which  compute 

necessary  tables  of  integrals , or  which  produce  tables 
of  points  describing  deformed  shells,  or  which  translate 
solutions  from  [6]  to  solutions  to  the  same  problem 
as  defined  in  this  paper. 

The  second  purpose  of  this  appendix  is  to  describe 
and  illustrate  SECANT.  This  description  is  especially 
non-rigorous  and  should  be  understood  as  a bag  of  tricks 
rather  than  as  scientific  fact. 

SECANT  was  modelled  on  Polak  [^3  by  Rauch,  imple- 
mented by  Marz,  and  modified  again  by  the  present  author 
in  Rauch  et  al  . It  is  designed  to  be  used  with 

a subprogram  GAUSS,  modelled  on  the  IBM  SSP  program  SIMQ 
[^3  » which  solves  systems  of  linear  equations,  and  with 
a subprogram  RESIDU  which  calculates  the  goodness  of  a 
guessed  solution  to  a system  of  equations  by  substituting 
the  guess  into  the  system,  and  taking  the  square  root 
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of  the  sum  of  the  squares  of  the  differences  between 
the  left  and  right  sides  of  each  equation  in  the  system. 
If  this  square  root  RESID  were  ever  zero,  the  guess 
would  be  a genuine  solution  to  the  system.  In  practice, 
it  is  usually  sufficient  for  RESID  to  be  less  than  an 
agreed-upon  small  number. 

The  subprograms  SECANT,  GAUSS,  and  RESIDU 
communicate  through  a block  of  COMMON  data  and  through 
argument  lists.  In  the  programs  below  the  author  has 
restricted  the  common  data  to  items  of  REAL  type  and 
the  argument  lists  to  items  of  INTEGER  type.  Except 
for  variations  in  the  COMMON  block  and  in  the  argument 
lists,  and  in  SECANT'S  DIMENSION  statement,  SECANT 
and  GAUSS  are  immutable  from  application  to  application, 
and  are  listed  below  only  once.  RESIDU,  on  the  other 
hand,  must  be  redesigned  for  each  application  to 
produce  the  differences  between  the  two  sides  of  each 
equation  (stored  in  one  column  of  a two-column  matrix 
E)  and  to  produce  the  residual  RESID. 

SECANT  starts  off  with  a central  guessed  solution 
and  tries  to  improve  that  guess  in  two  ways.  First  it 
takes  additional  guesses  by  modifying  eaoh  coordinate 
of  the  original  guess,  one  at  a time,  by  a fixed  amount 
EPSILO.  If  any  of  these  new  guesses  produoes  a smaller 
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residual,  that  new  guess  becomes  the  center  for  further 
guessing,  and  we  say  we  have  a local  variation  improve- 
ment. At  the  same  time  that  it  is  taking  these  guesses, 
however,  SECANT  is  building  up,  through  the  matrix  E, 
a table  HBAR  of  approximate  first  partial  derivatives 
for  use  in  creating  a better-aimed  guess  of  a Newton- 
method  type.  After  an  initial  period  of  relying  upon 
local  variation  alone  (while  the  table  of  approximate 
partials  is  being  developed),  SECANT  alternates  a 
local-variation  guess  with  a Newton  guess,  and  does 
not  accept  a local-variation  improvement  unless  the 
Newton  guess  of  the  same  iteration  fails  to  produce 
an  improvement. 

We  note  that  each  time  an  improvement  is  adopted, 
the  central  guess  is  changed  and  the  table  of  approx- 
imate partials  for  the  next  Newton  guess  will  have 
most  of  it 8 columns  thrown  slightly  out  of  alignment. 
This  potential  difficulty,  however,  is  usually  no 
difficulty  at  all,  and  indeed,  once  the  Newton  guessing 
starts  producing  improvements,  it  usually  continues 
producing  improvements  rapidly  until  the  central  guess 
is  so  good  that  the  old  EPSILO  is  too  large  to  produce 
meaningful  partials  anymore. 

The  key  to  efficient  use  of  SECANT  is  an  efficient 
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subprogram  RESIDU.  In  applications  of  any  complexity, 
machine  time  is  spent  chiefly  in  calculation  of  residu- 
als and  data  for  partials.  In  the  real  programs  which 


follow  the  program  illustrating  the  use  of  SECANT, 
every  non-obscuring  device  available  to  the  author 

: 

to  save  time  in  RESIDU  has  been  utilized. 

It  is  also  vital  that  RESIDU  not  be  called  unneces- 
sarily. That  is  the  reason  why  no  Newton  guesses  are 
tried  until  the  table  of  approximate  partials  has  been 
filled  with  reasonable  data,  and  why  the  Newton  guess 
is  only  tried  once  and  not  repeatedly  cut  down  and 
re-tried  as  Polak  suggests.  That  is  also  the  reason 
why  DTEST,  which  helps  control  the  size  of  EPSILO,  is 
kept  large,  and  why  ZTEST,  which  tests  the  size  of  the 
residual,  is  kept  flexible,  so  that  guesses  which  are 
good  enough  shall  be  accepted  as  solutions  promptly, 
and  the  machine  shall  not  waste  its  time  in  negligible 
and  expensive  improvements. 

The  real  programs  which  use  SECANT  have  built  into 
their  main  programs  and  their  RESIDU  subprograms  the 
ability  to  handle  continuations  of  both  types  described 
in  Section  VI  above;  this  feature  has  proved  economical. 

A feature  that  has  proved  wasteful  has  been  the 
ability  to  handle  both  odd  and  even  computations  (as 

J 
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defined  in  Section  VII  above).  Storage  requirements 
for  the  tables  necessary  for  odd  computation  are 
extravagant  when  even  computation  is  planned;  and 
features  in  RESIDU  that  save  time  for  odd  computation 
can  waste  time  for  even. 

We  proceed  to  the  example  used  in  the  illustrative 
program.  It  is  proposed  to  solve  the  system 

2 2 

x - pxy  + 5y  - 2xz  ■ 30 

2 

w - 2xz  + pwyz  - -20 

2 5 

pwx  - 3z  + zy  - - 19 

3w^  - 2xw  + 7z^  a 96 

In  thi6  system,  p is  a parameter  whose  initial  value 
is  2. 

This  system  was  created  artificially  around  the 
solution  w ■ 4,  x ■ 1,  y - -3,  z - 2,  but  we  start  off 
in  the  main  program  with  an  initial  guess  of  0,  1,  0, 

•5  • (In  the  main  program,  M is  the  number  of  equations 
and  variables,  ZPERM  is  a vector  that  holds  the  guess, 
and  PARAM  is  p.)  What  happens  during  the  run  is  illus- 
trative of  the  strengths  and  weaknesses  of  the  procedure. 

SECANT  fails  to  solve  this  initial  problem,  stopping 
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after  172  calls  to  RESIDU  at  w - 4,30,  x - 1.27, 
y - -1.13,  z • 1.93*  which  is  associated  with  the  large 
residual  33*24  . SECANT  then  passes  control  back  to 
the  main  program  which  does  a type  one  continuation 
on  the  parameter  PARAM.  For  PARAM  - 2.01,  and  with 
w,  x,  y,  and  z starting  at  4.30,  1.27,  -1.13,  and  1.93* 
SECANT  succeeds  in  31  calls  to  RESIDU  in  obtaining  the 
solution  w - 3*93,  x - 1.11,  y • -2.98,  z • 2.03,  which 
is  associated  with  the  small  residual  .0000890  • The 

continuation  now  proceeds  to  find  equally  good  solutions 
in  equally  few  calls  to  RESIDU,  for  as  many  new  values 
of  PARAM  as  time  permits. 

When  the  solution  PARAM  - 2.01,  w - 3*93,  x - 1.11, 
y - -2.98,  z - 2.03  is  inserted  as  the  new  starting 
point  in  a later  run,  and  the  continuation  is  directed 
backwards,  the  following  solution  is  obtained  for 
PARAM  - 2.00  — w - 3.96,  x - 1.10,  y - -2.98,  z - 2.02, 
residual  » .0000438  • Ue  are  surprised  to  find  this 

distinct  bona  fide  solution  in  such  close  proximity  to 
the  known  solution  w » 4,  x ■ 1,  y ■ -3,  z ■ 2 — the 
solution  we  had  hoped  to  find.  It  nay  be  that  some 
combination  of  continuations  of  first  and  second  type 
can  take  us  from  the  solution  we  have  to  the  one  we 
wished  to  obtain.  But  it  must  be  remarked  that  in  a 


less  artificial  situation  we  might  have  no  idea  that 
another  solution  was  nearby. 


I LL  US  T i<  A T l v H PKo  RHf.M  --  USE  OF  THE  SUriKQUIlNE  SECANl. 


C MAIN  ROUTINE 

C THIS  HOUT1NE  WEnO s IN  THL  INITIAL  DATA,  CALLS  SECANT, 

C AN.J  MAKES  ULCISlONS  AHOUT  WHETHER  AND  HOW  TO  CONTINUE 

C urjtt  S L C rti  I T SUCCEEDS  OH  FwlLS  IN  ITS  WORK. 

COMMON  ZTENP(IU) .ZPLWM(IO) .DEl  IA.WESIR, 

C E(2« 10 ) .FF ( lO.ll > .OTEST ./TEST.PAWAM 
C IF  THL  OIMLNSION  Of.  7PFKr1  IS  ( K ) » THEN  THt  DIMENSIONS 
L OF  / f Li"iF‘ » t,  AMt)  FF'  AHt  (M).  (2.M),  ANU  (M.MfI) 

C RESPECTIVELY. 

C 

C INPUT  UATA  HLWt. 

2Ptwr.(  l > = o 

7PLW*  ( > =1 
<PtKi*'(4>=ii 
JPLR.M**  > = .*> 

PAKAM=2 

C COMPUTE  OH  Ht AO  TN  SPECIAL  CONSTANTS  HEHE 
C (IHIS  PROBLEM  HLOUTRES  NO  SPECIAL  CONSTANTS.) 

C THE  MAIN  LOOP  BEGINS  here 
1 2)tST=.lt-4 
UTES  Tr.feF-OM 
LEL  T A = l 
DO  14  1=1. M 

14  IF  (OLLTa.LT .AHS<7PERM( 1 ) ) ) BELT A=ABS ( ZPEKM ( 1 > ) 

LALL  SECANT (M) 

CALL  l<EPO<<T(M) 

C A PARAMETER  IS  \)  a H T E 0 AND  THE  PROGRAM  IS  CYCLED  BACK. 
PAWA.rsPrtKAH-F.nl 
IF  (PAWAM.LT . (P.m  GO  TO  1 
STOP 
CM. 


I 


n n 


L. 

c 

c 

c 

c 

c 

c 

(. 

c 

c 

c 

c 

c 

c 

c 


ILL  Ub  Tit  A 1 i V/l.  PKOGK  fti'i  --  USE.  OF  THE  SUIEKUL1  TINE  SECANT. 


SUbHOD  f 1NL  S L C A M T < H ) 

This  routjnl  ATTEMPTS  in  solve  the  .system,  of  equations. 
COhMON  / T t MP  (10)  .?PLRM<  10)  .DELTA. RESID. 

C L ( 2 . 1 U ) » FF  (10.11  ) •DTE,  ST  • Z T EST  * PAR AM 
HEAL  Nil 

DIMENSION  OMEN A ( 1 0 ) , E TE MP ( 1 0 > . D ( 20 ) , DEL ( 1 0 ) . MB AH  f 1 0 , 10 

) .HTEST (10) 

the  dimensions  op  omega,  etemp,  and  del  ape  always  the 

SAME  AS  (HAT  QF  ZPtPH.  THE  DIMENSION  OF  L)  IS  TuICE 
1 H A I OF  /PLPM,  AND  THAT  OF  HHAR  IS  THAI  OF  /PERM 
SNU AhEd • THE  riiMt-NSION  OF  R TEST  IS  ALWAYS  10. 


C 

C 

c 

c 

c 

c 

c 

f 


MEANINGS  OF  VAnIAH|tS. 

/PL K m Is  THE.  VFCTOR  CONTAINING  THE  UEST  GUESS  SO  FAK 
F OK  THE  SOLUTION  OF  THE  EOLATIONS.  II  IS 
ASSOCIATED  WITH  THE  WES  I DUAL  WPERM. 

/TEMP  IS  THE  CURRENT  GUESS  VECTOW.  ITS  RESIDUAL  IS 
K TEMP . 

OhEGA  IS  A VECTOR  WHICH  CONTAINS  A BETTER  GUE^S  THAN 
T r it.  LUWKFNT  /E‘EWf1.  THIS  VECTOW  CAN  ONLY  HE  FOUND 
lit  A LUCAl  VARIATION.  OMfuA  IS  NCT  ADl  UMAT  jCALLY 
StW.STITUTEn  FOR  /PERM.  SUCH  A SUBSTITUTION  ONLY 
OCCURS  IF  nEW I ON • i METHOD  DOES  NOT  ALSO  PRODUCE 
IMPROVEMENT  IN  Tile  SAME  ITERATION.  THE  RtSIUUAL 
ASSOCIATED  WITH  OMEGA  IS  HRTEMp. 


t'  A Mi  ETEMP  CONTAIN  VALDES  OF  THE  LEF  7 -HANU-S I UES 
OE  TEif  M EQUATIONS  TO  UE  SOLVED.  E ( 1 . 0 ) IS 
ASSOCIATED  WITh  /PERM,  E(2.J)  WITH  ?TLMP«  AND 
ETEMP  WITH  OMEGA. 

DEL  AMU  HEAR  CONTAIN  FIRST  DIFFERENCES  FOR  USE  IN 
NEWTON'S  Me" THOU. 


C 

C 

C 

C 

C 

c 


FIRST  THE  PARAMETERS  OF  THE  COMHINE  D LOCAL- VAR  I AT  I ON 
-ANO-fvEW  TON-METHOD  ARF.  INITIALIZED.  SOME  OF  THIS 
IN IT  1 A Li/A  T IQI  IS  DONE  R Y THE  MAIN  PROGRAM  ALREADY. 
Tut  PARAMETERS  SO  PREPARED  ARE  DELTA.  /TEST 
and  DTtST. 

MOUNT  sfl 

WRITE  (t.,19)  DEl  TA.KODNT 
19  FORMAT  ( 7X  , 'DEL  r A IS  NOW  * . F10 . 3 • 6X . 1 5 > 

DO  21  1=1.10 
21  WTE ST < I >=1 .UEGU 
NU=lU 00000 
hSIJhsf) 

J=H  To 

IFLAG  = 0 \?>*V  F,,r 


1 


— — 


J 
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iLLLSTKAllVt  PROGRAM 


I'Sf  ti#-  THt  SUiiKGL  1 1 WE  St  C AN  T 


— 


; I 


1.0  2(>  1=1  .M 
1)1  1 > = ) 

C wl  HtGlN  hr  LSI  APL  r SUING  DATA  FOR  THt  FIRST  GUESS, 
ion  do  no  iri.ii 
110  /I  t IP  I X >=ZPEKMI  ) 

CALL  RESllHJI  ltK) 

KOIJ.4  T sKOI.IL  T ♦ 1 
FPERMrKt  SIO 

WK1IE  |fc.l'4U4>  PPF  KM.KOUN! 

IF  CHPEHiS.LI  .ZlfST  ) 00  TO  1400 
C THt  I TER A T I ON  HEGInS  HERE.  FIRST  COMES  AN  ATTEMPT  TO 
C IMPROVE  THL  CURRENT  GUESS  HY  VARIATION  OF  A SINGLE 

C VAR1AHL  F.  IF  THIS  SUCCLtl.'S « THF  IMPROVEMENT  IS 

C SI  OH ED  TEMPORARILY  IN  OMEGA  AND  ITS  ASSOCIATES. 

C A r I HE  SAME  TlYt.  UATA  IS  PE  I AG  GATHERED  FOR  THF 

C dENEFIT  OF  NEWTON'S  METHOD. 

200  IF  <J.EG.2*M)  J=0 
J=J*1 

AOO  tPSlLiJ=AMNl  IDE  1.1  A.  Ml) 

400  JJ=J-  (U-1»/y|*H 
DO  4 OS  I z 1 .11 
4 Oh  l IENP  < I ) = ZPLl?M  ( T 1 

/It  MP4  JU)=7PLRMl  JJ > ♦FPS 1 LO*f)  < J) 

CALL  RESIGN!  2.M ) 

KOUN  T =KOL'NT ♦1 

RTErtHrutsn 

IF  IHTEMP.GE.7TFST ) GO  TO  hOO 
ZPt'RMI  JJ  )-7.  TEMPI  JJ) 

HPEKMsP 1 1 MP 
GO  TO  140U 
500  DO  hlO  lsl.ii 

DEL I I ) = <K I 2. I I -FIT . 1 ) 1/EPSILO 
510  h<MMl  I.UJIsOEL  I T >*U(  J) 

HOO  IF  (RTFiIP.GE.RPERM)  GO  TO  6*0 
IFLAG=1 

t>2u  DO  630  1 = 1. M 

L TEMPI  I )=£I2.T  ) 

630  OMIGAI 1 TsZTErtPI i ) 

RKTEMPsRTtMP 

b<#0  IF  IHOUWr.Lt.fi)  GO  TO  1200 
C NOW  NEWTON'S  MCTHOD  IS  GIVEN  A CHANCE  TO  TRY  TO  IMPROVE 

C I HE  GUESS.  IF  IT  SUCCEEDS.  IIS  SUCCESS  IS  GIVEN 

C PHECEOENCF.  EVLN  OVER  A LARGER  SUCCESS  OF  THE  METHOD 

C UF  local  variation,  if  it  fails,  the  LOCAL>VARTATION 

C IMPROVEMENT  SIOHFU  in  OMEGA  --  IF  THERE  was  INCHED 

C A LCCAL-VAHIATTON  IMPROVEMENT  --  IS  ADOPTED.  ELSE 

C NO  change  takes  place  in  zperm.  in  any  of  these 
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ILLUSTRATIVE  PKOGHAM 


■HHHHEEET 


USt  GF  I ML  SUItHULT  lNt  SECANT. 


C THREE  LasES.  rut  1IEKAIICN  rit&lNS  AGAIN  AT  L I Nfc  200. 

7 0U  00  /In  Isl.ii 
10  71n  1 1 = l . h 
710  FF  ( I , l I ) = h 1 1 A H ( I , 1 1 ) 

00  720  1 = 1 «M 
720  Ft-  ( 1 » M ♦ 1 )=F  ( 1 . D 

CALL  GAUSS(M.KS) 

IF  ( K S • F.  (j  . 1 I Go  TO  119A 

40o  do  yin  1 = 1,11 

910  2 TF.hh  I J >=*PEKK  J >-FF  ( I , M + l ) 

CALL  RESIUJ(2,M) 

K0UNT  = K0lifgT  + l 
KTF.  MP=hES10 

IF  (NIK  rip  • G T . ? T L ST  ) GO  TO  1000 
DU  “3  20  1 = 1.11 
420  2PEKM ( 1 ) -l  TEMP ( T > 

rplhm=k  rt  mp 

GO  TO  1 A 0 0 

1000  IF  (RTtfIP.Gr.RPFKM)  GO  TO  1200 
1FLAG=0 
GO  1010  1=1  ,H 
/.PERL  I I ) -L  rt-’iP  ( T > 

1010  E ( 1 , I ) =E ( 2 , I ) 
rpehm=h  tl.ip 

GKIlt  (fa,  1403)  RPt  KM  » KOUIV  T 
MO=0 

DO  102  0 l = 1.ri 
1020  NU=iJU+FF(  1 ,M  + 1 >**2 
f\JU  = NU  * * • K 
bO  TO  200 
119R  *SUn  = rfSJMl 

whIIc.  (G.1199)  kSU.1 
IF  ( A SDH. tO. M I GO  TO  1400 
1194  FORMAT  (/X.'bAOsS  H A 1 L E f J 0Ur*,6X,Ib) 
l*rOu  IF  (J.LI.2+M)  GO  TO  1300 
DEL  T A = DtL  T A/2 

IF  (CEL rA.LT.DTFST ) GO  TO  1400 
wKIlt  <o,19)  DELTA, KOUNT 
1300  IF  ( IFLAU.F.d.O)  GO  10  200 

1 F L A G = 0 

00  1310  I = 1 . M 
E ( 1 « l ) =E  TEMP  ( T ) 

1310  /HF.KM(  I )=OMEGA(  T ) 

KPEKMsKKlLMP 

Wi<  1 1 E (fa,  1402)  WPEKD.KOUNT 
DO  1330  1=1,4 

1330  P TEST  ( 1 L-l  )=KTEST(  10-1  ) 'ippV  ~ " " 

RTEST  ( 1 ) =KPtRM 


illustrative  program 


list  OF  THE  SUBROUTINE  SECANT 


IF  <H  ft  iT  (U/PTFSTi  10  ) .U  . I ,«#9)  ) GO  TO  200 
C THFkE  AWL  SLVTKAl  WAYS  FOK  lhl.S  SUPKOUTTNF  TO  TERMINATE. 

C A SUCCESSFUL  TERMINATION  OCCUKS  IF  RPEHN  GETS  SMALLER 

C THAN  ZTLSF.  UNSUCCESSFUL  TERMINATIONS  OCCUR  UHLN  DELTA 

C GETS  TOU  SMALL.  OR  WHEN  TEN  LOCAL  VARIATIONS  IN  A HOW 

C PKUUUCE  LESS  TFAN  A 1 PER  CENT  CHAN&F  IN  THE  RESTOUAL. 

C OR  WHEN  GAUSSIAN  REOUCTION  FAILS  TO  PRCOUCE  A TEST 

C VECFOR  FOR  NEW  TOMS  METHOD  FOR  THE  M-Th  TIME. 

14(10  Will  IF.  lo,l4U4)  RPEKM.KuUNT 

1h02  FORMAT  (fcX.E14.7.FX, IS.6X, ‘LOCAL  VARIATION** 

1403  FORMAT  ( GX . L14 . 7 .FX , I 5 . 6X • • NEW1  ON • * 

1404  FORMAT  ( toA  . 1 14 . 7 .f>X  . 1 f>  * 

RE  TURN 

END 
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IlUJSIKAI  i VE  (■’  l<  l J t>  t<  A I • 


DSL  OF  THL  SUHKOUT  INE  SECAM 


SUHROul  ll\L  GAUSS  ( D « K S > 

c this  routine  solves  ststems  of  linear  filiations, 
common  i ri.MHd  u > ./permi io>  .del ta.klsio. 

C C.CP.l'j)  .PH  I II.  in  .RTF  ST. /TEST  • PAR  AM 
C GAUSSIAN  ME DUCT  I ON.  WITH  PIVOTING. 

MMrrtO 
KS  = U 

T UL= . iE-1 b 

C FORWARU  SOLUTION.  WORK  On  THE  J-TH  COLUMN. 

00  6b  0 = 1. M 
Jj=J*l 

H 1 G A = 0 

C FIND  THE  ROW  WHOSE  >J-TH  COEFFICIENT  Is  LARGEST  • 

DO  30  I=J.M 

IF  (AMS<HxGA).Ar4S<FF(  I.J)))  20.30.40 
20  PIGA=FFU.J) 

1 iT  A A — I 
in  CONTINUE 

C IF  THE  DIGGLST  J-TH  COEFFICIENT  IS  TOO  SMALL. 

C THtN  G 1 OK  UP. 

IF  (AHSU.lGA)-lol.)  3b.  St).  4 0 
■b  rS=) 

HE  TURN 

C MOVE  I HL  ROW  WITH  THL  BIGGEST  J-TH  COEFFICIENT  TO  THE 
C U-Tn  POSITION.  REPLACE  IT  RY  THE  FORMER  J-TH  ROW. 

4 1)  00  Hb  K = I .'lil 
SAVL'sFFI  J.K) 

FF { J,K )=FF( IMAX.K) 

Hb  FF( IMAA.K)=SAVL 

C UIVIUL  THL  J-TH  ROW  UT  ITS  LEADING  COEFFICIENT. 

UO  SO  *\  = 1.MM 
SlJ  FFl  J.K)=FFI  J.K  >/BlGA 

C ELIMINATE  THE  J-Th  VARIABLE  FROM  THE  RF.MA INING  EQUATIONS 
IF  (j.tN.M)  GO  TO  65 
00  HO  i=JJ.h 
00  60  KSJJ.HM 

60  FF I I »K  ) =FF I I .K  )-FF( I ,J**FF (J.K) 

6b  CONTINUE 

c hALKWANO  SOLUTION.  COMPUTE  THE  N-TH  VALUE  INTO  FFlN.M+l) 

00  70  J=2.M 
JJ=M. 1*J 
JJJsJ-1 

(10  70  I = i.JJJ 

1 I=1M-I 

70  FF ( JJ .M •« ) =FF ( JJ ,MM ) -F  F ( JJ . 1 1 ) *FF ( 1 1 .MM ) 

METURN 

END 
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11  LUS  I wM  ] v/t  1‘KOliHAPi 


USt  OF  ThL  SUHkOLTINt  StCANl. 


suuHouriML  hlsiduiname  .m> 

C IHJS  HCltrjixK  F'KtPUHEK  A KfSlUOAL  AWL*  1 HF  AKKAY  E F OH 
C THL  SOitHOUl  IMF.  SFCAMT. 

co. amok  l rt>'P<  lu  » »7pehm<  10 ) .uely  a .kesio. 

C LI^UD  .PF  (10*11  ) .OTEST  ./1EST  .PAH AM 
\»-i  TEMPI  1 ) 

X = <?T£  MH  ( t ) 

Y-ZTcMP ( A > 

/ztTl'.t*V  ( A ) 

E FIMAME  . 1 ) = <*t?-pAKAM*X*t  + 4»r**2-i»*X*?-30 

eimmf:  )=*-«!*  x^*f,ahaf,*w*y*z**2*<'u 

t (NAME.  3)=PANAn«-b*X-3*?*2*Y**2-X**t>-19 
t «UAM[.  «H  )=3*'***:>-;?*X*W*7*/** 3-9f 

rtCSlO=E  < M'rtfAC  .1  » **?n  ( NAME  «2)**2*E(NAME  , 3 ) • + ?*F  { M AME  . H ) 

**? 

HtSU>=hESIU**.b 

hctohn 

CNO 


■nwiryssa—si 


1LLUSTPA1  I \it  PKOCjKAM  --  OSL  OF  THL  SUHPOLUNt  SLCAmT. 


SUtjhCiU  T 1 ImL  UtPOMKM) 

C THIS  HOuTIMt  WrtlTtS  HP  1 Hi.  OUIPUT. 

COMiil'N  2 T L MP  < 1 0 ) • ?PC  PM ( 10)  .UtLTA.MLSIU. 
C t (2»1»)  .FF«J  n,  11 ) ,0I>  ST  .71EST  tPAKAM 
VKITF  I o)  in  ) p/w<AM 
10  FUHi'iAT  < // . 10X  « »fAHAM  =*,E14.7) 

00  2 0 1 =1 «M 

20  ifKin:  (b./^n  t./piki-md 

21  FOHI1AT  (10X.i?,^)(,Ll4.7) 

via  ft  (o,22» 

22  FOKMAf  ( // . lx ) 

PETUHfj 

tlMLi 


NAJUrt  PrOGRAM 


null  UuLkLIWo  OF  AN  LLLlFSOiuAL  SMLLL 


C 


C 


C 


C 


MAIN  HOOT  I ML 

COMMON  DELTA, »Lt  , bLLSN  « A l TCH  , j<HO  • GNU  t ZTES1 • DTE  ST  , 
C CHAO  *A0  ) »C?<  AU«A() ) ,(.  A i 40  I .SGi A0«A0  ) ,S7<  TO  ) e 
C sue  A i),  A Ut  AO  ) .901  Ag,  AO  ) ,E<  2.401  «FF  ( AO  «A1 ) . 

c z rtrviuo)  ./perm  Am  ,puci  .o<ao> .HPLRr .resiu 

READ  IN  DATA. 

REttO  lb.*/)  XX.CKEHEN 
IoT  rt=x X 
HEAU  <b.27)  X 
27  FORMAT  <ElA.7.febX> 

M = X 

HEAD  <b.Z7)  HH0,GMU.PEE»AI TCH 

Z ILS1  = . 7bt-l)A 

OTtoTs.ME-OA 

DLlTa=  . 0001 

JFL AG  = 0 

STAHIsliLLlA 

HEAL)  Ib.ifM)  (H(T)*1=1,M) 

?.t  Format  <a( tx.fia.7) » 

1 HK[SU  = I«L) **2 
CALL  HEAOlN(n) 

MAIN  LOOP  STAR  IS  MFKF  . 

2 DO  AO  1=1 .M 
AO  ZPF KM( 1 IrHI 1 » 

IF  (lOTit.LT.l)  GO  TO  A1 

T H'iPOKsRMu 
HHO=ZPtRM  IOTA  > 

ZPLHM < IOT A)=1EMpOR 
A 1 CALL  SFCANTIM, IOTA) 

IF  I lOTA.LT.l)  00  TO  M2 

TtMPOMshllU 


A 2 


13 

A 

211 

?3h 


RHO=ZPERM(  mi  A ) 

ZPEHr  i I01A)=TFMP(JP 
IF  (HHFRf  .G|  , I .MF.-U*))  STOP 
CALL  ORAPMIiI.JFlAG.IOTA) 
JFLAGsl 

ZTESl=Af1AX7  I .‘SvrPERM.  .7bE-UM  ) 
OLLTArSI ART 


UO  7 A 1 = 1.0 

IF  ( AHSfP ( I ) ).LT .1 .1F-A0) ) HIT  ) =0 

PARAMETER  TS  WAHlEU  AND  THE  PHOGRArt  IS  CYCLED  BACK. 
IF  (IOTm.LT.  1)  GO  TO  238 

HI  10  I A)=HI  IOTA  HCPtMEN 

00  TO  2 

khosrhoachemln 

GO  TO  i 

tNU  copy  mm*  to  rr r 


5* 


yhu-v 


n n n n 


MaOOK  --  f>Ul>  HULxLlNG  OK  AN  ELLIPSOIDAL  SHELL. 


SUUKUUl INL  KtAUlMM) 

niMrtUK  UE  L T A « Ht  F « HLE.SU  * A I l C H .WHO « GNU < 2 TEN T .UTCST  . 
C Cl(4(lt«iM.C?(un»40)«rA(4  0».S#i<‘‘U»AO),S7(40>. 

C StU  *0  .40  .4U  ) « $4140.41)  ) »E  (2*4()  ) iFK  140.41  ) . 

C /TEMP (40  ) .2PF  RM(40)  .P(40>  .H(  40  > .KPEHM.RESIU 
THIS  ROUTINE  READS  IN  THE  INTEGRALS. 

Till  LuOP  L l M T T S AHL  SPECIAL  TO  THL  CASE 
M=40  IN  TMl  PROGRAM  ••  INTEGRALS"  WHICH 
PRODUCES  THE  VALUE*  BEING  HEAD  IN  HERE. 

DO  200  1 = 1.40 
HEAD  (<).2Ul>  T ,r3I  I ) ,S7»  I ) 

FORMAT  ( 12.E14.7.E14. 7.50X) 

00  210  1=1.40 
00  210  J= 1.40 

REAO  (b.20?)  T . J.C1T I .0) .L2T I « J) «Sbl I .w) «S9< 1 . J) 
FORMAT  t 12. 12, El  4.7. E 14,7. E 14 . 7 . E 14 . 7 . 20X > 

00  220  N = 1 ,Pt>7U 


200 

201 


210 

202 


READ 

(«. 

203)  11 

c 

14 

. J4 

, K 4 . V 4 

21S 

s«m 

. Ji 

. t\  1 > = 0 1 

Stt(  1 1 

. K 1 

,0.1  >=01 

Sn(dl 

.11 

• K t ) = 0 1 

S«<01 

» K l 

. 11 >=W1 

SfUKl 

. I 1 

.01 >=0l 

So  ( K 1 

.01 

, ll)=\/i 

21b 

SB  I 1 2 

.02 

,r\2  > = 02 

SO  (12 

«*2 

,02  >=02 

S6(02 

. 12 

, t\2  ) = 02 

S*(  J2 

. K2 

, 12  )=V2 

Sb(*2 

. 12 

. 02 ) = V? 

SH(K2 

. «j2 

. 12>=02 

217 

S8(  13 

.03 

. a 3 ) =03 

Sb  (13 

«K3 

• 03 ) =03 

Stt(03 

.13 

, K3 )=03 

saio.i 

• K 3 

. 13 >=03 

SbTK3 

.13 

, O 3 ) =03 

Sb(K3 

.03 

.13) =03 

?lb 

Sb(  14 

.04 

»K4 )=04 

S«(  14 

.K4 

.04  )=V4 

Sb(  J4 

. 14 

, K 4 ) =04 

SOT  JH 

.K  4 

. 14  >=04 

S M ( K 4 

.14 

. 04 ) =04 

Sii  T K 4 

.04 

. I 4 ) = V4 

220 

CON  I 1NUL 

203 


FORMAT (4( 12. 12. T2.L14.7I  I 

RETURN 

ENU 


jn/Mi  *r>i  r ?n 
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D; 


MAOUK  PHOOKAr.  --  oul)  DUCKLING  OF  AN  tLLllSOlLAL  SHELL. 


SUoKOUl  l UN  RlSillUINNN.fi.  IOTA) 

COMMON  .)»  LTA.Mtf  .HtLSO, Al ICH.RHO.NNU.ZIESI  .01FST, 

C Cl  I'Ki^Ol  •C?i4U»40)  »C3 ( 4 0 ) »SU (40,40)  » S7 ( 4 0 ) « 

C SO(40»*»r»»4U>  «SS(40.40>  *k  {2i40>*FF(40iH1)  • 

C ZTtNP<40  ) ,2Pt  RM40  ) ,P(4U  I ,0(40  ) »RPF HP  . RES  10 

c this  routine  calcoi ates  hperm  and  eii.j)  cr  it 

C CALCULATES  KTEMP  AMO  E<2.J>  FOR  THE  BENEFIT  OF  SFCANT 

C THE  VARIABLE  K«k  TURNS  OUT  fj  BE  1 IF  Ifj  IS  O0n. 

C.  2 IF  1*J  iS  EVEN.  WHEN  THE  DEEPEST  LOCP  STARTS  FROM 

C KKK  Al'O  SUES  BY  THUS,  MANY  USELESS  NULT IHLICATIONS  BY 

C ZERO  ARE  ELIMINATED,  AND  THE  SPEED  OF  THF  SUBROUTINE 

C IS  DOUBLED. 

IF  ( IoTa.LT  .1  > c;0  10  41 

Ttf’iPORsKHD 

kmu;/ TEMP ( IOTA ) 

L ft  BP  ( 10  1 A ! = Tri“toCR 
4\  CALL  PrOBPlM, IOTA) 

FAC  TCRS12*  ( i-r,N|j»*2  ) *BEESO/  ( AITCH**3 ) 

un  iio  lsi.n 

tPHEflSO.I) 

IKAUSsO.il 
DO  10  4 J=1  ,r) 

TEBPsy .0 

K*K=2-I  14J-U  I+ki)/2)*2) 

no  ior>  n=akn,m»? 

IDS  TEMPs  rtNI’«SM(  f , J,K ) *P(K  > 

THAUSslhAuS^Cl ( T«J)*H(J) 

10B  EPMERstPHKH-SM 1 , J) *P< J)-< ,5*HH0*SS< I * J ) -TEMP ) *0 ( J ) 
110  El  NUN, I )=-TKAMS*F ACTOR* <-.5*RM0*S7 ( I > ♦EPHEM) 

KESIUSO.U 
DO  120  Isl.Pl 

120  RtSinsKESIO+ElNNN.I >**2 
R t S 1 OsRt  S 10 * * . S 
IF  (iulA.CT.l)  HtTUKN 
TEdPORsKHO 
RH0=21EHPI IOTA  » 

2 temp i iota)=tfmpor 

Kt  TU*U 
ENi> 


nprv  IV  Ml  t T° 

ubfl  r'  p • • 
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MAJOR  PROuKArl  --  <)OU  tiuCKLING  OF  AN  ELLIPSOIDAL  SHELL. 


SUHKGUTlNL  PLOMPIM.JOTA) 

COMMON  01  LT  4 »HLF  .BLESO. AITCH.KHO.GNU.ZTEST.UTEST, 

C CU4Q*4  0I.C?<4U.4n)*C3(40>«Sb<40'40ltS7(40>« 

C S<1  <40.40.40  ) »SV(  40 .40  ) .t  I?  *40  ) »FF  ( *M)  »41  ) t 
C i f t>iP  ( 4 0 ) * ZPE  Mh ( 4 0 ) » P I 4 0 ) » H ( 4 0 ) .HPERtf.RESIU 
C THIS  ROUTINE  COMPUTES  THL  P»S,  GIVEN  THE  B‘S.  FOR  THE 
C BLNtFII  OF  THF  ROUTINES  RESIOU  AND  GRAPH.  FOP 

C THE  MEANING  AND  I'SE  OF  AKK,  SEE  THE  COMMENT  FOR 

C THe  ROUTINE  PEsIOU. 

00  b J=1.H 

5 nil ) =7  tempi  I > 

NO  lu  1 = 1 .'1 
GO  10  J = 1.M 
10  FF(I.JI=C«iII«J) 

00  20  1=1,0 
20  FF  I 1 .M-U  )=RHO*CM  I > 

00  bU  1 = 1. M 
EPHEMsO.il 
DO  2 7 J = 1,M 
TEHPsft.il 

KKK  = 2-I  I-fJ-I  I T*.l>/2)*2> 

00  2b  KSRKK.M.2 
2b  TEMPsTt  rtpMStJI  I . j.K  »*HIK  > 

?!  E P HE MsEPHEM* I Sb { I.J»-.b*TFMPI*8l Jl 
30  FF(l.K*n=FF(T.M*l  > ♦ A I TCH*HEESG*tPHE  M 
CALE  GANSSIM.KS) 

00  40  1=1. 0 
40  PII»=FMI.rt+l) 

RETURN 

ENO 


M/iJOH  PNiJtiKi  I -I 


<)00  OllCKLlNG  OK  AN  ELLIPSOIDAL  SHtLl 


SUUHOUT IwE  GRAPH (Mt JF  LAG. IOTA  ) 

COMMON  UtLT  A «HEp .HtLSO. A I TCH  « RHO  » GNO . Z TEST  « UTLS»T  • 

c ci(4n.uij)tc?(uo.4n),r3<4  0),s6i‘»Otio>tS7i‘*o). 

C SHI40.4U.4U)  tS^(40»40)  • E ( J?  » 40  ) » FF  (40*41  ) * 

C /TEMM40)  ,ZPEhK(40)  .P(4U).U(40>  .NPERM.RLSID 
DIMENSION  XPOS(Si!)  .YPOS(F2l 
C KtPOHT  (Ht  SOLUTION  A NO  SOME  OF  THL  PARAMETERS. 

WHITE  (uill)  OlF  , AlTCH.GNU.  IOTA 
11  FORMAT  (//tlX. • H E E = *,E14.7**«  A I TCH  s ••El4.7«*« 

GNU  = *. 

C E14. 7 . • . IOTA  = • . 13) 
wH  1 IE  (0,14)  WPFHM 
14  FORMAT  ( I X,  'KFSTOUAL  =»,E14.7) 

00  lb  1 = ) .M 
ZTfc  iilJ(  1 >=l«t.RM(  I ) 
lb  HI  I Is/HtKM  I ) 

CALL  HCOMPIH. TOTA ) 

WHITE  (b.lb)  (H(l)tl=l,M) 
lb  FORMAT  <4< IX.E14.7)  ) 

C CALCUIaTL  THE  new  POSITIONS  OF  f'ift  y-thkee  points  on  the 
C ELL  IPSO 10.  IF  JFLaGsO.  which  is  so  the  first  time  THIS 
C HOUT INF.  IS  CALLED.  REPORT  THEM.  IN  ANY  CASE  HEPORT 
C RHO,  POL  AH  DEFLECTION,  AND  EQUATORIAL  CONTRACTION. 
IFLAbsO 

FAC TOR  = A 1 1CH**3/(12.0*( 1.0-GNIJ**2)  ) 

HATI0=1  /HLTSQ-1 
Fsl./AI TCH 
L=s 

PI=3. 1415*3 
R*=PI/(b2*L> 

LL  = bi>L 


1000  w0A  = 0 

no  ion  j=i,ll 

XI=J*0* 

CO=LOS( XI  ) 

S I =s In ( a I ) 

OEf\lOfcsHEESO*CO*CO*S  I*SI 

SOKTUsSOHl (UENOM) 

XsSl/SWRTO 

Y = -l)EESU*CO/SOHlO 

H0ALPH=SI4HATI0*SI**3 

K2  = S1 ♦Sl/UENOM 

BE1 A=O.U 

PS 1=0.0 

HilPsO.O 

DO  10  I=1.M 

CS=COS(  UXI  ) 

SWsSIM  1«aI  ) 


n*\PY 


dim’  iwr 


■fn 


r* 
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MAJOR  I'KOtKAfl  --  ODD  HuLKL  1 Mj  OK  AN  LLL  IPSOIOAL  SMt.Lt  * 

I 

ML  TArHE  I A ♦ H ( 1 >*sN 
PSI=PS1*P(  I )*SN 
PSIP  = PSIPM*P(  1 )»CS 

io  continue 

PHI  = X 1 -Hi  r A 
SP  = SII'J(PHl  ) 

CP=C0S ( Pm 1 > 

K2lKCW  = -.S*HHO*Sl*i>I/Uf  NOK 

1 EMP=SP/ (KOALPH«AlTCH) • (X*AI JCH+PSI *CP*RZEKOV*SP-KOALP 

H*GN(j*PSIP- 

C K2*r,NU*HHO*SP)-BEFSy*Sl/(SQRTU**3) 

TLrtP=IEIlP*OA 

WOAsWOrtFTtMP 

IP  (J.bT.  ( J/l.  )*l  ) 00  TO  100 

XPOS  ( J/L  >=X>F*  (qOALPH*PSIP-R2*RHU*SP-GMj'»PSI*CP-GNU*KZ 

EROV*  SP ) 

YPOS  t J/L ) sY  + wnA- TtHP/2 
ion  CONTINUE 

IF  tIFLAG.EN.lt  C-0  TO  108 
LtFLt  C = t.FL  - YP0Stb2> 

CONlHA=l.-XPOS(?fc) 

WHITt  to, 104)  KHO.OEFLEC, COMMA 
104  FORMAT  t 1 X » • MHO  ='*t14.7«*  INFLECTION  =*.F8.5«*  CONT 

HACT10N  =• 

C « F tl  « b ) 

if  t jflag.ne.u*  pe town 

HP  I ft'  (0,107)  yPOS(52) 

107  FOKMAT  ( / . 1 X , • UF  FORMED  SHAPE!  NORTH  POLE  AT  • 

« F8 • 4 ) 

108  IF  (lFLAb.EN.lt  WRITE  (8.104) 

104  FORMAT  (lX.'llFTAlL  NEAR  SOUTH  POLL*) 

00  110  1=1.17 

Jin  WRITE  (b.Ul)  XU0S(  I ) .YPOSm  .XPOS(  17*1  ) ,YPUS<  17-H  ) .XP 

OS(34+l>. 

C YP0S(344I) 

111  FORMAT  (1X<5(F*).4«2X.F8.4*4X)  ) 

WRITE  (8.200) 

200  FORMAT  (//.20X) 

OX=PJ/I270**»L> 

IFLAb=lFLAG>l 

IF  (1FLAG.LT. 2)  GO  TO  1000 

ML  TURN 

END 
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OAOOH  PKOGiIaK  --  f VLN  bUCKl  1 NG  OF  AN  LLLlPSOiUAI.  SHELL. 


C Mn  1 N UOUTIML 

COiMlOf,  DLL  T A • MtF  »PLI.«s(J  , A 1 1 LH  .MHO  »CNU  • / TES  T » OTLST  . 

C Cll  20.20 > .L2(20.2l) ) .C3(2<>  I .S6(?U»?0I  .S7I20)  , 

C S#(20.2U.20)  '7(2  0.20  ) .E  (2»?0  ) .FF  (20,21 1 . 

C l rEPP^O  > .2PLHM20  ) , P ( 20  I « H ( 20  ) .MpEKP  .HLSIO 
C MEAD  IN  UATA. 

PEAU  (bt?/l  XX.r»<EMCM 
lUf  A=<x 
KEAU  (5.27  ) X 
27  FOWMAT  IFJi4.7.6AX) 
i'l  = X 

NEAi)  (b.2/»  MHO.GNU.DLE. AITCH 

2IESI  = .7f>L-0H 

orESTi.bE-nH 

DLL T A= , OUUl 

jFLAftrf, 

SI ArtlzJELl A 

PEAL)  (b.?U)  (R(T).l=l.N> 

2tt  FOMilAT  (HI  IX.F1H.7)  ) 

1 nEESQ=I«tt*»2 
CALL  HLAOIAKM 

C MAIM  LUMP  STARTS  HfKE. 

2 DO  HO  i - 1 • M 
HU  2PLKr1( I l=M( I I 

IF  (iOIA.LT. II  CO  TO  HI 

ILf'iPOMrPhO 

HKO  = XPLH^( TOT  A ) 

ZPLKM  IOTA  ) = TFMPOK 
HI  CALL  SFCANT(M.IOTA) 

IF  (I0TA.LT.1I  GO  TO  H2 

TLMPOP=KMO 
KH0s2PLH0( IOTA  I 
ZHEWf' ( iOTAlsTFdPOH 
H2  IF  ( HPt'PO.G  T . ( . SL-O?  ) ) STOP 
CALL  GH APH ( (1 , JF  L AG  . IOT  A I 
JFL AGsl 

2!EST  = ArtAXl(  ,5*BPEH(1«.7bE-0H  I 

OELTAsSTAHT 

DO  13  i=l.rt 

7.X  IF  ( Ar»S(H(  III.LT.I  .IE-301  t H(I»sO 
C A PAHAdLTF.H  IS  VAIIlED  ANO  I HE  PROGRAM  IS  CYCLED  PACK. 
IF  (iCfa.LT. II  GO  TO  23# 

237  HUOTAlsm  IOTA ) ♦CKLrtEN 
GO  TO  2 

23#  HH0sHHO>CHE#EN 
GO  TO  2 

END  .,>3^ 
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o o n n 


SUUKOU  I INL  KLAuTN(F') 

tUtl.-lUl'l  DELTA. HtFI  .UEESO.AI  TCH  . l<HO  » GNU  . 2TE  ST  , (JTF  ST  , 

C Cl  120,20)  ,02(20,20)  * C 3 ( ? 0 ) , S 6(20*2  0)  , S 7 < 2 0 ) , 

C SB  <2 1)  ,20,20)  .sy<  20.20  > .1  (2«20>*FE(20.21>. 

C /1ENP(2U>  .ZPERtM2n  >.P<20),B<20)  . KPERf’ . HE'S  ID 
C THIS  PUUTJuE  Ht*nc  IM  THE’  INTEGRALS  . OhlTlING 
ALL  THOSE  WHICH  Have  ANT  CUD  SUBSCRIPTS. 

THL  LOOP  LIMITS  ABE  SPECIAL  TO  THE  CASE  tf  = 40 

in  the  pkoghan  "integrals"  which  pkoolces  the 

VALUES  BEING  KF  At)  IN  HERE  . 

1)0  200  1 = 1, *+0 
HEAD  (B.2D1)  J.X.XX 
IF  ( (J/2)*2.LT.J>  GO  TO  200 
J = .J/2 
C3( J)=X 
S / < J ) =AX 

200  LUUriNUL 

201  FOHrtM  <T2.tl4.7»El4.7»b0X) 

DO  21o  1=1,40 
DO  210  J=l,40 

READ  (0,202)  li.Jj.X.XX.XAX.XXXX 
IF  ( ( J I /2  ) *2,1.  T.  II  .OH  . ( JO/2  ) •2.LT.0J  ) GO  TO  210 
11=11/2 
JJ=JJ/2 
Cl ( II.JJ)=X 
C2( II ,OJ)=XX 
S6( I l . Jo)=XXX 
SS( II , JJ ) =xxxx 
210  CONTlNUr. 

202  FOkrtAT  ( 1 2 . 1 2 . c 1 4 . 7 , £ 14 . 7 . E 1 *♦  . 7 . E 1 4 . 7 . 2 0 X ) 

00  220  N=l.2ri7U 

HEAD  (4,203)  Tl.Jl.Kl , V 1 , 1 2 . J2 ,K2 . V2 . I 3 . J3 , K3 . V3 , 

C I4.J4.K4.V4 

IF  ( ( 11/2  ) *2.L  I , 11. OH. ( J1/2I*2.LT. Jl.OH. (Kl/2 ) *2.LT ,K1 

) GO  TO  216 

11=11/2 
J1=J1 /2 
Kl=Kl/2 

5H(II . Jl , X 1 >=01 
SHI  11 ,M.Jl)=Vl 
S«(  Jl.IlfKllsVl 
SB( J1.K1.I1 )=V1 
SB ( K 1 « 1 1 . J 1 >=V1 
StHiM,ol,H)=Vl 

216  IF  <(I2/2>*2.LT.12«0R.(J2/2)*2.LT,J2.0H.(K2/2)*2.LT»K2 

1 GO  TO  217 


N 

| 

r 


12=12/2 

J2=J2/2 
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MaJOH  HUobMAM 


tVCU  HULKL  ING  OK  AN  ELL  IPSO  I UAL  SHLLL. 

I 

K2=*\?/2 

S8(I?.d<f,K?»sUi» 

S U I irr.  ,Hd,U2)  = \/  2 
SUI  J2,12,K2)=\/2 
S8 I J2 » K2 » I 2>=V/2 
su(K2.i2.j  2>=v2 
SdO\2.  Jd,  1 2)=\/2 

217  IF  (<13/?»*2.l.l.l3.0K.IJ3/2)*2.Lr.J3,0H.lK3/2>*2.LT.K3 

) GO  10  216 

13=13/2 
j3=J3/2 
K 3 = K 3/Jr 

SUt I3*J3.K3)=V3 
Sfc  I I3.Ki,J3  )=\/3 
Sdl J3. 13,H3)=V3 
Sd(o3»K3»l3)=V3 
SUIK3. 13, J3»=V3 
Sfc(K3,J3,i  3)=V/3 

21 6 IF  ( ( 1*4/2  U2.LT  . I**. OH.  ( J4/2U2.L1  . J4  .OR.  (K4/2  >*2.LT.K4 

) GO  TO  220 

J4=I‘»/2 
J4=JH/2 
K4=K4/2 

SU<  1*4  « J4.HH  ) =WH 
Sftl  IS  ,<<*♦,  J4  >=\/H 
Sfa (JS , IS.KS )=VS 
su i jh , i< s , m»=vs 

S«(K'4,IH,o4)=V4 
SU(K4,J4, I4)=VS 
220  CON T 1 NIJL 


203  FOWrtMCM  12.12, T2.tl4.7U 
KE1UKN 
LNU 


MA  JON  PKUOK  All 


FVtN  BUCKLING  Of-  AN  ELLIPSOIUAL  ShLLL 


I 

I 


! 


SUuNGUT  I lit  NESIPUINimN.M.  101  A) 

COMMON  DEL  TA.BEF  .HtE.SO  . A I TCH  » KHO  • GNU  . ZTE  ST  .UTFSl  . 

C Cl(20.2U».C?<:>0.20l.C3(?0>«S6<20.pn>,S7(20J, 

C Sb(?0*?0.20  > *s9(2u*?0 ) «E  <?«?n ) «FK (?0«?J ) • 

C 2Tt>IP<2U  ) .7PF«M<20  » ,P<?  01.6(20)  .RPERN  .RES  ID 
C THIS  ROUTINE  CALCULATES  ((PERM  ANO  t(l.J)  ON  IT 
C CALCULATES  KTFNp  AND  E(2.J>  FOR  THE  HEAEFIT  OF  SECANT. 

IF  (10Tm.LT. 1)  pO  TO  <41 
TEMPONrKHO 
KHU  = Z TEMPI  TOTA  ) 
l r E.MP  ( IOTA  »=TFHpOR 
41  Call  PCOMPJM, IOTA) 

FACTCRsl?*(l-fiNlJ**ii)*HEESO/(  AITCH**3) 

do  no  (=i.ri 

EPHEMsO.O 
TNANS=0.0 
00  106  Jr 1 ,h 
TEMPrO.O 
00  105  Krl.M 

105  TEMP  = TLi'1P4S6(  T . J.K  ) *P(K  ) 

TNANSrTHttNS  + Ll  ( I » J ) *R  < J ) 

106  EPHtNrEP*lt*l-S6<  I . J)*P(  J > - I .5»HH0*SS < I . J ) -TEMP ) *B ( J ) 

110  Ed'INN.  1 )r-TKAMS*FACTOK*(-.5*NHO*S7(I)  ♦EPNEM  ) 

HESIOrO.O 
DO  120  1 = 1. M 

120  KLSIO  = HLSIO>E ( NNN. I T * *2 
HESIU=HESID**.S 
IF  (IOTA.LT.l)  HETUMim 
T£MPON=NhO 
KHOs-ntrtPI  TOTA  ) 

2 TLMP ( IOTA)=TFMPOR 

RE  TURN 

END 


[;;)py  jiyijj  ,r«p  Tp r ■ • : 
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o n 


MAJOR  PKObKrth  --  EVLN  RUCKLING  OF  AM  ELLIPSOIDAL  SHELL. 


SUnKOUl  i ML  PC(>.-IP<M,  IOTA  ) 

COMMON  ME  L T A.FlL  •HlLSU.  A I TCH . RHO . GNU . ZT EST « OTfST  , 

C Clip  0.20 ) *C?  < ?0»?0  ) ,C  3<  20  > . S6<  20 .20  > .S7  (20), 

C S8120.20.2U)  .54(20.20)  « L ( 2 « 20 ) « FF  ( 20  < 2 1 ) • 

C /TEhP(2U>  .ZPEpM<  2(1  ) ,P<20)  ,H(20>  .RPEKN  .RES  ID 
THIS  KOUriuj  COMPUTES  ThE  P*S.  GIVEN  THE  H*S.  FOR  THE 
OENEFIT  OF  THE  ROUTINFS  RfcSIUU  AND  GRAPH. 

oo  b r=i ,m 

b U < I ) =2  TE  MP  < I ) 

r.iu  m i — i * m 

00  10  J = 1.M 
10  FH  I . J)=C2<  l . J) 

uo  20  I = i . n 

20  FF(I»M*l)=MMn*C.MT> 

00  40  1 =l.m 
LPHLNrO.il 
UU  27  J=1,M 
r E K P = 0 . L* 

00  25  Krj.rt 

2 5 TEl-.P  = TEMP  ♦SOI  I . J.K  ) «K(k  ) 

27  FPHEF=EPHEM-MS6< 1 » J ) - . b*TE KP ) *B ( J ) 

30  FF ( I ,M+l )=FF ( I » M ♦ 1 ) ♦ A I TCH*«EESG*LPHEM 
CALL  GAUSSdl.KS  » 

00  4 0 1 = 1 . M 
40  P ( I ) =FF  ( I « fi  + 1 ) 

RETURN 

LNO 


' ■'V 
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MAJOR  PROGRAM 


EVEN  DUCKLING  OF  AN  ELLIPSOIDAL  SHELL 


SUbHOUT  1NE  GRAPhTM. JFLAG, IOTA ) 

COMMON  UElTA.HE  E.BEEStt.Al TCH ,HHO .GNU . /TEST . DTEST . 

C Cl(?n.?(j)«C?(?n.?0)«C3(2ft)*S6(?Q«?0).S7(20)t 
C SB <20. 20.20)  .s'*!  20,20  >.F<2.?0  >.FF(PO,  21)  . 

C /TEMP  <2U  > . /PE  OM ( 20 ) .P(20).H<?0)  .RPERM.RESIU 
DIMENSION  XP0STK2) . YPOS ( 52 ) 

C REPOKT  THE  SOLUTION  ANu  SOME  OF  THE  PARAMETERS. 

WRITE  » MEt.AlTLH.GNU.  IOTA 

11  FORMAT  (//.IX.'hLE  = *.E14.7»,»  AITLh  a *»E14.7**» 

GNU  a * . 

L E 1 4 . V » • . IOTA  s*.IJ) 

WRITE  (A. 14)  RPTRM 

14  FORMAT  ( lx,  »HFSTUUAL  =*,E14.7) 

UO  IS  1=1.4 

/tempt i )=zperm<  t ) 

15  H(I):m.RhU) 

CALL  PCOMHTM.TOTA) 

WRITE  (fe.ln)  (M(I).Irl.M) 

1 G FORMAT  (4T1X.E14.7)  ) 

C CALCULATE  THE  NEW  POSITIONS  OF  FIFTY-THREE  POINTS  ON  THE 
C ELLIPS010.  IF  JPLAt=U,  WHICH  IS  SO  THE  FIHST  TIME  THIS 
C ROUTINE  IS  CALLED,  REPORT  THEM.  IN  ANY  CASE  REPORT 
C HHO,  polar  ijeflfctton,  ano  FOUATORIAL  contraction. 
IFLAGsO 

FACTORS A I TCH**3/(12.0*| 1.0 -GNU** 21 ) 

RAT  10  = 1 / Hf  FStt-1 
F=1 ./AITCH 
L = ? 

PI=3. 1415*3 
DX=PI/(52*L> 

LL=52«L 
1 ()  0 U W 0 A = 0 

UO  100  Jzl.LL 
XI=J»UX 
COsCOST  XI ) 

SI=SIN(X1 I 

UEN0h=HLEStt*C0*C0lSI*SI 
SORTDsSnRT  (DF.NOM) 

X=SI/SnRTU 

r=-BtFStt*C0/Si4RTU 

R0ALF’H=SI4RATIU*SI*»3 

R2=S  I *SI/OENCM 

PETASO.O 

P S I = 0 . 0 

HSIPsO.O 

UO  10  1=1. M 

CSsCOST I*XI ) 

SNsSINT 1*1(1 
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MAJOR  PROGRAM 


FVtl'J  RUCKLING  OF  AN  ELLIPSOIUAL  ShLLL 


Gt  I a = HL  I A + H ( i ) *«;N 
PSJ=IJSI+F<  J >*S|V 
l,SIH=HSiP  + ,-**l*M(  11+CS 
10  CONMNUL 

PFI 1 = XI -OF.  T A 
S P=Glw(PHl  ) 

CP=C0S (PHI) 

RZF.KO  V = “.  Ei  + KFHT  + J?  1 *S  I /(>F  Nt)M 

I LMP=SP/  ( KOALPH+Al  I CM  ) * ( X*A  l TCH  + PS I * CP  + RZEROV * SP-HO ALP 

H*Uf\IL*PSIP- 

C H2*GNU*KHu*<:p)-BLt;S(l*SI/(SNKTD**4) 

TFIMP=  IfMP  + nx 
UOA=W(lA  + TLM» 

IF  < J.G  I . < J/L  ) *(  > GO  TO  100 

*POS< J/L )=<+F* ( WOALPH*PS1P-H2*RHO*SP-GNU*PSI*CP-GNU*RZ 

tROV*SF  > 

YPOS< J/L >=Y+W0A-TEHP/2 
100  CUNT  I NOE. 

IF  ( lFLAG.Fo.l ) GO  TO  lOfl 
OEFLFC  = H*-t  - yPOS<52) 

CONTRAsl ,-XPOS ( ?6 ) 

WHITE.”  ( G « 1 0 4 I KMO.OLFLf  C,  CONTRA 
104  FORMAT  < 1 At  1 RMO  =',tl4.7,*  DFFLLC  T I ON  =*.F8.b,'  CONT 

K ACT  ION  =• 

C • F a . b ) 

IF  (jFl.AG.l'iL.in  RETURN 
WRITF.  (b.in/)  yPOS  ( 52  ) 

10  7 FORMAT  ( / « 1 X t * OFFOKMt  D SHAPE!  NORTH  POLE  AT  * 

* F B • 4 1 

100  IF  (IFLAG.LO.il  WRITF  (6,109) 

109  FORMAT  (lX,*l’ETAlL  NEAR  SOUTH  POLE*) 

00  110  1 = 1 ,17 

110  WKITt  (G.llll  *pOS( 1 ) , YPOS(I) ,XPOS( 17+1  ) , YPOS( 17  + 1 ) ,XP 

OS( 34+1 ) , 

C YP0S(  34  + 1 1 

111  FORMA  1 (1 X,4(Ffl.4,2x,F8.4.4X)  ) 

WRITE  (b,^00) 

200  FORMAT  <//,20X> 

DX=PI/(27U4  + l.) 

IFLAG=1FLAG+1 

IF  (1FLAG.LT. P)  GO  TO  1000 

RETURN 

ENU 


— 


hAJOK  PKOGBrN  --  oUl>  blJCKLlNb  OF  A SHhLH  1 CAL  SHELL . 


C MAIN  HUH  I 1 Alt 

COMMON  H|iUS.AOM«ONU.KOOl  • I J f. L TA.ZltST.Hfl  KN.Kl.SIM. 
c HLOE  . L LA  OH  .0  TF  ST  .1  ABLE  ( Ah • 35 • 35 1 . E ( 2 . AS ) « FF  < 35 . 36  » « 
C ZTE-IR  ( Ah)  .ZPEkM  Abl  .AI3SI  .HI  3b  > 

C HLAU  IN  DATA 

HE AO  (b»2/>  X.CrEMEN 
IOTA=X 

HEAD  <S.2/>  X 
Ms  X 

HE  aO  <S.p7)  RMOS.GNU.AOH 

OELTAs.OOGI 

KOIJM  f = 0 

START  =ntl TA 

GTEST=.bE«OH 

ZTESTs.PbE-OH 

HEAIJ  O.?0>  (A(I)«1  = 1,M) 

H001  = ia.-i.*GNl!  **?)**.  5 

KCUE  = <?.  / ( AoH*AUH*KUOT » ) • ( 1 . ♦GNU/ I 2 . * AOH*ROOT 1 ) 
FLA0H=1  . ♦GNU/  ( 2 . »AO||*KOOT  ) 

?7  FORMAT  (f 17.10.E.SX) 

2b  F OHM  A f (4<  lX.FlH. 7.5X1) 

C REAM  IN  TABLE 

00  JO  I=l.Sb 
00  10  J=1 .3b 
oo  io  k si, a*> 

10  TABLE ( I . J.KJsO.n 
HEW  I NO  rt 
00  20  L=1 . 1 0650 
REAM  (0.21)  J.J.K.VALUr 

IF  (I.GT.35.0R.J.GT.35.0H.K.GT.35)  GO  TO  20 

TABLE I I • J.K  JsVALUE 

TABLE  I l *K.J>sWA(  UF 

TABLE  Id. I «K  JsVA|  UE 

TABLE  Id. K • I I sVAl.UC 

TABLE (K.I.J)sValOF 

TABLEIK.J.I  )=VA|  UE 

20  CONTINUE 
ENU  FILE  ti 

21  FORMAT  II2.IP.I9.E14. 7) 

C THE  MAIN  LOOM  STARTS  HEKE'.  SET  UR  ZPlHM  AM)  CALL  SECANT 
1 00  HO  Is]  , N 
Hu  ZPEKM  1 IsAli  I 

IF  (iOTA.LT.ll  GO  TO  SO 
TF  MPsRMOS 
HHOSsZPEHrtl IOTA  ) 

ZPEKM I IoTa  JslFMp 
50  CALL  SECANTCM.IoT* > 

' IF  I IOTA.LT.l)  GO  To  60 


j’flpY 

UUf  I H 3 •‘.if  M 


• i.iti'id.  it — , inu^wuilwii 


i'iAJUK  PHUOMrtri 


nun  HllLKL  IM]  UF  A SPHLH1CAL  SHLLL 


TLMHsKHOS 

KH0:>  = 2PLHF<  IOTA  ) 

2PFKK( 10TA»=lf rip 
60  CALL  OUTPUT  (M.IOTA) 

C ML  “SL  T SLC  Aim  T *6  PAAAMtTLKS 
KUUuTsO 

IF  iMPLM^.GT. ( ,?bi-03) ) STOP 
2HST  = anA*II ,5*MPEKm.,P5E-0H ) 

CLL T ArS I AMT 

C A PAMAWFTLK  IS  VAMlEU  AMO  THF  PHOGMAh  IS  CTCLCU  HACK. 
IF  I 10TA.LT.lt  GO  lu  200 
A(IOTA)=A(lOTA)*CKLrLN 

GO  TO  1 

20U  HhuSsprtUS+.OO* 

00  TO  i 

1 NO 


Ujr  • 


4 S i 


V - K 
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HAJOK  KHllLKAri 


OUU  UUCaLING  OF  A SKHtKICAL  SHLLL 


SUbBOD  I i Nt  COMFM  (M.IOTA) 

COBrtON  KMUS.AOH,  GNU,  HOOT,  DLL  l A.ZTLST  , HPiKM . KES 1 0 , 

C KCOF » LL  A OH « n TrSiT « T AHLF-  (35t35«3b)  .F  < ?. . 35 ) • FF  ( 35 , 36  I • 
C 2TtHP<  3S)  tZPEuM  AS) . A<  3b>  .ti<  3b> 

C COMKllTt  T Hk.  H*S.  THt  LOOK  LIMITS  AKL  St  T UK 
C JO  AVOID  UMWFCJ.SSAKY  MULTIPLICATIONS  HY  2EKO. 
KAT10=A0M*B00T/(  LA  OH  ) 

DO  2 l = l.li 
^ aii  ) = /rFmM  1 1 
DO  «dOu  M=i  »n 
TlVl=?*lv>l 
I >VI'il  = 2*N*  ( M + l ) 

T = U 

DO  10H  Lrl.ii 
IA=I AKS(l-N) 

IH=l.b*.b*J-l)**IA 
IC=BAXO( I A « IB ) 

IOsMlNO (M.L+N) 

DO  inO  KsICtlD,? 

101)  T = I ♦ TAHLC(L,K.iV)«TNl*A(L)«A(K)/TNNl 
?0  0 B(N)=BAT10/IN*(M-H)-l-GNU)*JT/?-A(N)  ) 

KLTOKN 

END 


ggpy  n"u; 

^rPTFir  rr  . 


M(\JUK  PiUHjHmM  --  ODD  lilJtKLlNli  UK  A SPHLRICmL  SHELL. 


SOhkUUT  1 Nt  KuSlOU  ( MMf-i  . M « 1 01  A ) 

COMMON  HIlOS  . AOh « ONu » ROOT  .DEL  T A • Zl  EST  , RPE  K«v«  . RES  10  , 

C K L Of." « t L 0 OH  * D T f S T »TAHLt43b«35».5b)  « E 4 2 * 3 5 t . E F ( 3 5 « 3 b ) • 
C L I E -'.P  4 3b>  .ZPEqXifi)  . A4  lb)  .H(  3*j> 

COMPUTE  THE  RESIDUAL  AND  THE  VECTOR  E FOR 
THE  USE  OP  IHt  S4JHR0UT1NE  St  CANT.  THE  LOOP 
LIMITS  AKt  SET  UP  SO  THAT  UNNECESSARY 
MULTIPL ICAI 10HS  HY  ZERO  MAY  HE  AVOIDED. 

IK  (IOTM.LT.lt  DO  TO  l 

TEMPsilHUS 

F<HOS  = ZTEMP(  IOTA  ) 

Z I E M P 4 IOTA  )=TE,“IP 
1 CALL  COMpb(M.IOTA) 
f-0UK  = 4 ♦rti>uT*tL  Al)h*AUH 
no  2 0 0 Mri.M 
TNl=2*N*l 
TNNl  = 2*i'J*  ( N+l ) 

T = 0 

DO  lUU  Lsl.M 
I A = 1 A H S 4 l.-hl ) 

IB=l.b*.b*  4-1  ) **IA 
IC  = MA<i|4  I A » I M ) 

IOsMlNO  4 N.L  + N ) 
bo  100  KslC.IIl.? 

100  T = 1 ♦ 1 AMLE4L,n.N)*TNl*A<L )*P4K)/TNN1 
if  0 0 E 4 NNN  . N I = 4 4N*  4 N*1  ) - 1 + GNU  ) /F  OUR -MHOS  ) 

C * l\  4 0 ) _2*B  ( N ) ♦?• T 

RES  1 U = il • H 
DO  300  1 = 1.  M 

.-SOU  HES  I0  = KES I O+E  ( Nf\jN . I )**'/. 

HtStl’  = l<ti>lH**.i> 

IE  4 1 0 1 A . LT.it  RETURN 

TcMP=RHOS 

RHOS  = Z TEMPI  IOTA  t 

ZTENP4  IOTmIsTF.HP 

RETURN 

END 


MAJOR  PHOGKAU 


--  nun  BUCKLING  OF  A SPHERICAL  SHELL. 


SOBBOUT 1NL  OUTPUT  (M.IOTA) 

COflrtUN  KMOS,AOH,GNU.POOT.[)KLTA.ZltSr.WFEHM.KLSIll. 

C KCUF .f  L AOU.nirST ,T ABLE < 35 • 35 . 3b ) t E < 2 , 35 > » FF  < 55 . 5h ) « 

C *rtr,F<  Ab)  ,7P£RM(iS>  *A(  351  tH(3S) 

no  10  I - i 

A ( I »=/PLPh(  I ) 

10  2TErtP<  l »=2Pt.KM<  i » 

CALL  COt'iprt(h,  IOTA) 

WtUJt  (o.lt)  RHOSt  AOH.GfJU 

11  FORMAT  tlX.’HHOS  =*tE14.7.»  AOH  =».FA.3t'  GNU  =*.F 

6.  A) 

WIUTE  (c^ll  (A(1>»I=1,M) 

21  FORMAT  (4( lX.ri4.7t5X )) 

CALL  UEFtLC  <B> 

WL  tUKN 
END 


•Miyy  ivjm  imp . r ^ • f 
uur  * h • > 
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r/VJijK  *-*l< U»-»< A M --  nOi'  HOLM.  INI.  of  a sphlhical  sihll. 


SUIlROU  T 1NL  UFPLFC  ( M I 

tOi'i.'iOf  HUOS, AOH.bNU.HOOT .DLL  TA.ZTLST  . KPFHrt  . HLS ID , 

C KLUt  .Ll.AOH.nTF ST  . f AHLF  ( 4h  . 45 . 25  ) . F ( 2 . 4b  > . FF  ( 35 . 36  ) . 

C 2 T t KP  < 3b  > , 2PF RM  Af>  > . A < 3b  • . B ( 35  > 

C COMPlirt  POLAR  OF  Ki.LC  T ION. 

W=2-HHUS*KCUt * AOH* < .8*GN0*A(2)+(1-GNU) ) 

00  10  L=2.:wl,2 

111  W = W-2*<rt(L)-KC0F*AUH*U<L)*<  1 ♦GNU  » - C l -GNU  » *KHOs*RCOL*AO 

H# A ( L )/2> 

00  20  usj.H 

211  W = W-N*(0+l.)«A(W)**?/(?.*N4l.) 

T = II 

1)0  4 0 * = l.f1 

KrtlNUS=  ( l\-1.  ) - < (K-1  >/2»*2 
00  30  0 = 1.0 

A M 1 I\’IJS  = ( H”1  > “ ( (^-J  >/2)*2 
MT0P=(  U-l+N.'lirjUS  )/2 
KT  0P= (K“1+KM1PUR)/2 
U=-K*0*PH  1 Nl\l  ( K , N I 
DO  2 '}  JM.'iKII’ 

?i\  u=0-8  * ( H*  J-2*Art  ilvil'S-1  (iPKINMK  ,2*J“2*NF  IOUS-1  ) 

1)0  2/  0 = 1.8  100 

2 / U = J-Ai*  I 4*  J-2**MTNUS-1  ) *PK  1 NN  < N . 2*  J-2«KMMUS-  1 • 

00  2fr  i=1.KroP 
00  2L  J = 1 . 1’j TOP 

2p  U=U-  ('♦*  I-2»Krtl  f||iS-l  » *(H*J-2*MH1NUS-1  ) * 

C PH  lMi-i  I 2*  1 - ? •kM  NUS  - 1 .2*J-2*NPUNUS-1) 

3 IJ  T = T>U*U  < K » *rt(M) 

W = U + ( 1+GNU) *AOH*KCOt*T 
T = 0 

DO  <»U  K = 1.M 

DO  HU  A'  = 1 « ,v) 

NM1UIIS=  (M-l  >/2>*2 

N T i)P=  ( N-i  ♦ • IM  I Nu^  ) /?■ 

U=-P«INN(K,N) 

DO  3d  1 = 1 .AjTOP 

48  D=IJ-<H*i-2*NPfi\H|S-l)*PHlfMN(K.2*l-NMIMUS-l> 

40  T=T>U*MlK)*A<M»*K*<K+l)*N 
W=W-GWO*AOM*HCGF *T 
1=0 

00  bU  K=l,M 
00  b J 0=1.  M 

U=-PH  INrtUK-l.NM  •♦PH  I N M f K - 1 .N-l ) +PH  1 MN  < K ♦ 1 . h|4  1 ) -RK1NN< 

H*1  ,N-1 ) 

MPAHT  Y=0-<N/2>*3 
l\irOP=«IOI\)PAKTY)/2 
00  '*b  1 = 1. 0 I OP 

*♦«  U=1)^K*1K  + 1 )*«H*t-?*NPAm  r-1  »•(-?«  INN  < K-1.  2*  I-NP  ART  Y-l  ) 
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I'l  A JOK  PHOOKArt  --  ODD  huCKLINO  Of-  A SHhtH  I CAL  SHILL. 


L ♦FK [NNfK*l  ,2*l-HFArUY-l  > ) 

b()  T=  r-U*A  (k  ) *A»N)  *K*N*  (K+l  ) »<NM  >/(  <2*K-»1.0)»(2*N+1 .0  ) ) 

u=w-mihos*aoh*«coe*  r 

T=U 

DO  bU  (m  = 1 • M 

60  T=  !•*■?.  *t'J*fN*l.  I*AIN)*H!N)/|?.*N  + 1.) 

W=W*RCQE*AOH*T 

W / - - 'ft 

w i - w 2 / «i 

WHITE  (b./O)  M.uW.Wi 

7n  FOKriM  ( 1 a , •UTSTANCC  Of.  TUFF  A'  POLES  = ',F£.5. 

C /.IX, ‘TOTAL  POLAK  DEFLECTION  = ‘.FA.5. 

C /.lX,»SFfT  POLAK  DEFLECTION  = *,F6.b» 

WHITE  (b.OO) 

80  FOKriA  I f J of  oHX  ,#*•*♦,/)  ) 

HE  1 Urti'l 
LNO 


rnpy  mil  '-nr  rn  rr-r  f 
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hAjOK  HRObHAtt 


n'JU  HuCKLINb  OF  A Sf’HLHICAL  SHtLL 


I 


f UlILT  (ON  HIUUW(K«N> 

c INNtK  HH«.UIICI  OF  L^Gk-NUKt.  FUNCTIONS  K ANU  N. 
FKjNNsO.il 

IF  (IV.Nt.Kl  RFHlHN 
PR INw=2 . / ( 2.*N-M  . ) 

Rt  TURN 
tNO 
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MAJOR  BROOK)*  --  t.  \Jl  N BUCKLING  OF  A SOUtRlLAL  SMELL. 


' 

C HA  IN  K()tJ  1 I wt. 

C Oil  MON  KUOS  , AOH  , GNU  , OOO  T , BEL  Trt  . / TEST  , KPE  Kri  , Kt  S I P • 

C t <2 .3»  ) .FF  < 3(1  . XI  ) .?TLMP(  3ll ) . /PERM < 30 > ,A(30),R(3n), 

C I ABlE  130.30.30)  .OIEST.RCOE.ELAOH 
C READ  IN  UATa 

HE  AO  (S.2M  X.CREMLN 
10  f A = X 

READ  < rj  « 2 / > X 
M = X 

HE  All  KUOS.  GNU,  AOH 

DEL  TAz.UoUl 
IJ  TEST=  .R0t  -OH 
START=OElT A 

7Tti7=.^f>t-0«r 
READ  <M.?0)  ( A ( T ) » I=1.M> 

HOOTr  ( 3. -5.* GNU* *2  > *♦  .5 

KCOt  = (2.  / ( AOH*AOH*KOOT ) ) * (l.*GNU/<2.*AOH*ROOT) ) 

ELAOHr  1 .-.GNU/  I 2 . *AOH*ROOT  ) 

27  FORMA  r (t  l7.lll.G3X  ) 

2 A FORMAT  ( 4 ( 1X.F14.7.SX ) ) 

C READ  IN  lAliLf 
11  DO  10  1=1 .0 

00  Id  J=1.M 
(JO  10  k =1.1-1 

10  TABLE  ( I ,J,«  )=(!.() 

REMIND  <1 
CO  20  L=l,10«50 
HE AO  (4.21)  I .J.K, VALUE 

IF  ( 1 ,GT.2*M.0U.J.GT  .2*M.0M.K.GT.2*h)  GO  TO  20 
IF  ( ( I/2)*2.LT.I.0H.(J/2)*2.LT.J.0K.(K/2)*2.LT.K)  GO  T 

0 20 

1 ABLE ( 1/2. J/2.K/2)=V ALOE 
I ABLE  ( J / <r.K/2,J/2)=V/ALUE 
TAhl E (J/2. 1/2. K/?)=V ALOE 
TABLE (J/2.K/2. 1/2)= VALUE 
T ABLE (K/2. 1/2. J/2)=VALUE 
T ABLE (4/2. J/2. 1/2) =VALLE 

20  CONTINUE 
EKU  FILE.  4 

21  FUMMAT  ( 12. 12,  I2.E14 . 7) 

C THE  MAIN  1 000  STARTS  HERE.  SET  UP  ZPf.HM  AND  CALL  SECANT, 
l I'O  ho  i = i.n 

411  /BERN  ( I ) = A ( 1 ) 

IF  (lOtA.LT.l)  GU  TO  47 
TEMHslIHOS 
HMOS=/HtRM< IOTA ) 

2PEHM( IOTA )=TEMp 
4 t CALL  SECANT (M.lnTA ) . 
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M A J 0 K A I' i 


FVCU  HUCKl.  If.iG  OF  A SPHfcKICAL  SHklL 


IF  < I (J 1 A . L ( • 1 ) fiC  T(>  4 0 
T tnPsw hUS 
KHOSs^Pt PM (IOTA) 

( 1 1 ) T A »=Tf>p 

4v  call  ou thu r < m , iot a ) 

C Kt-Str  bfcC  ttM*S  PApAft  T Li<S 

IF  lKPt.Hh.Gr.  ( -11.1)  ) STOP 

ZTfc  S TrAlAXl  ( .5*RFfct<l%  .PSL-04  ) 

Uk  L T A = S T AKT 

C A PAKrtFit  TtH  IS  VAHTLU  AND  I hL  PKOGKAM  IS  CYCLtU  HACk. 
IF  I IOTA. I T.l  ) r,C  TO  SI 
A(IOTA)=MIOTA)  ♦ C.KLMEN 

bu  To  1 

SI  HHUS  = HHOS  + CI<tT't.fi 
go  ro  i 

k.  NO 


I 


NAJUK  PHObrtAil 


FVtfJ  bUCKL  1 UG  OF  A SPHERICAL  SHELL 


SUbKOU  I 1 WE  COPPn  (M.IOTtt) 

LOM.lCr:  KHOS.AOH.GNU.HOuT.OELTA.ZlEST.KPEKtf.RESIl), 

C J0>*Fr(5O«^l)«?-(LPH(3U)  ,2PFKM<30)  ,A(30).H<30)« 

C T ABLE ( Ml  « 30  « 30 ) • OTEST . KCOt . EL AOH 
C COMPUTE  THE  H'S.  1 HL  LOOP  L1MTS  ARE  SET 
C UP  TO  AvuiO  IJIJNF.LFSSAKY  MULTIPLICATIONS  BY  JEKU. 

HAT  10=AUH*RU0T  / ( 2*ELAOH  ) 

UO  2 l =1.M 
2 A(  I »=/TthP(  l ) 

00  2 00  i\l=l.H 
MIM=N-H\I 
T N 1 = 2 * N rJ  ♦ L 
TNNl=2*WN»  (NU  + 1 ) 

T = 0 

UO  10  0 L = 1 ill 
rA=lAHS<L“IVJ 
lC=FTAxo<  I A • i ) 

10=MU  0 ( h .L-Mvi  > 

UO  100  KriC.IO 

100  T = 1 ♦ TABLE  IL  ,K»N)  * TM1*A  (L  ) *MK  1 /TNN1 
200  b ( I!  >=RA  f lu/  (I\IW*  ( IMN  + 1 ) -1-bKU  ) * ( T/2-A ( N > ) 

HE  rUHM 
END 


i 

1 


-ncy  I’-  * 


L 
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MAJOk  HHtMjiu.n 


EVEN  IUICKI  llJ(i  OK  n SKhtHlLAL  St  ILL  L 


SUuKLU  I 1 !>,L  Kt  S l OU  ( UNAi  «N  . 10T  A I 

LOM'iOM  khOS  , AOI,  ,GMU  ,0  00  T .nfL  T A .Zl  EST  ,HPE  |<M  .Kfc  $ ID  , 

C (•  { Z • 3 0 I »hE|iOiU  1 » / T E MK ( 3 0 ) «ZKtKM(30)  . A(3U)  .0(30)  . 
C 1 A Hi f (Ml . 30  » 30  ) .DILST.KL OL.LCAOH 
IK  (1014.1.  r.  II  00  TO  t>0 
C COMPUTE  T tiL  HtESIUljAL  AMU  THE  VECTOK  L FOK 

0 THL  iJSL  OF  THE  KKnGP  AM  SECANT,  I HE  LOOK 

C LIMITS  AKE  SF  I UK  TO  ttV/010  UNfot CE SS Art Y 
C .MULTIPLICATIONS  liy  ZLPO. 

TCciHsKIiOs 
HHOSrZ  r t-f.P  ( luTrt» 
i ( E <viK  ( 1 01  w ) = 1 F i*.p 
50  CALL  COrtpBTM.ICTA) 

FUlirtm  » rt  0<*  r * E.L  A()H*  AOM 

00  ZOO  i\i=l.n 
i>.i\l  = M*M 
TNI  =2*  Ain*  1 
Ti  tl'd  l = 2*U(\i*  ( .JH+I  ) 

T = o 

LO  10  0 Lrl.fl 

1 A=1 AMS  ( L -•  I > 

1C=MA*0(  ( A . 1 ) 

1 U=N INu  ( M .1  ♦U  ) 

LJ  100  K = iC.lD 

10  0 T r I + (A'tLf.  (L  ,K.IJ  ) *Tul*A(L  1*0  (K  l/TNM 
ZOO  t(i':NM»!.l  = ( <:\|N*  (Mf.-Ll  ) -HOMi ) /FOUK-NHOS  ) 

C * A ( I.  ) -Z*H  (M  r 

HLsin=o.o 

UO  300  1 = 1,  (A 

300  HESID  = KLS11'  + E (I\inN,  1 ) **2 
KESI J=0LSin**.5 
If  (IOIa.LT.I)  UtTUUU 
TfJlPsMT'US 
MHOS"/  T Ei'il'  ( IOTA  I 
ZTLHP( 1 01  A 1 =TEhp 
He  TUKN 
LWu 
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bUbROO  r 1 IgL  OUTPUT  (M.IUlA) 

COMMON  KUUS,  AOH.bNU.ROOT  ,UEL  TA.ZTEST  .KFtRM.MLSlP. 

C KI2*3i))«FM3n,^lt  . l l'Ff.P  < 30  » .ZPFRM ( 30  ),A(30).B(3n). 

C TaMLL ( 3u . 31).  30 ) .DTEST .RCOE .ELAOH 
DU  10  i = l . h 
A ( I )=/l»U<h(  1 > 

10  ZTE.AP(  Da^PLKfM  I ) 
call  comph <m, iota > 

WRITE  (b.ll)  MHOS  t AUM  » ItNlJ 

11  FORMAT  (IX, 'KHOS  =*»L14.7,*  AOH  s'.FtJ.S.*  GNll.  s*,F 

6.3) 

WRITE  (b.2l)  ( a ( I ) • I = 1 • ft ) 

21  FORMAT  (mi.X»Kl4.7ibim 
CALL  Or.FLEC  («> 

Ml  I UR '4 


MAJOK  HKOumm 


t'VtU  hi.JC KL  l WG  OK  A bPHt.KiC.AC  bHLLL 


SUhKUU  f IPt  OIK-lFC  (h) 

COWMOV  KUOS  « AOh  • itfaU  .MOOT  • OtL  l A . I TfS  T , KFLMh  • KES  1 1) , 
t t (2 . 3 J ) .E  F ( 3l>  • 3l  ) ./lt>P(  30  ) . ZPEKM  (3n>,A(30).8(3n). 

C TAULt:  < 311. 30 .30  » .UTEST.KCOt  .LLAOH 

I'mtlgf.u  s.stlp 

COT'.PU  T t POlAH  uFKiECriOlM. 

s=rt 

bit p=2 

w=2-Kuns*KC0t  * A OH*  ( . ft*UfJU*  A < 2/2 ) ♦< 1-GNL  > > 

UU  1U  L=2«S.2 

10  4;  = V -2*  { A(C/2  > -WrOr«AOH*B  I L/2  ) * ( 1+GNU  > - ( 1 -GNU  > *KHOb*  KC  0 

E*AOH 

C * A ( ( / 9 ) / 2 ) 

no  20  j-s rpp , s . s i f p 
2 0 (,|  = U -hi*  ( U+l  . ) * A ( (\i/2  ) **2  / ( ? . *N  + 1 . ) 
f = 0 

UU  30  K=s ICP.S.S IEP 
KHINUS=<K-1  >-(  (K-l  )/?)*2 
[’0  30  f.  = STFP.S.sTFP 
NMlNUS=<f'-l  )-(  ( M - 1 > / 2 ) * 2 
M70P=  ( hj-l+Wrtih'US  )/2 
*aPP=i*-\*KHU>rO/2 

U--K*iM«PK  Ihlhl  ( K . M ) 

no  2«  j=i.urop 

pH  u='J-K*  < A*0-2*MF  jNUS-1  ) *PM  TNN  ( K • 2*  J-2  *NM  NUS- 1 ) 
no  27  o=i.«rop 

? / U=U-A  » ( J-2*Ki*;tMUS-1  ) ♦ PH  INN  (hi  ,2*  J-2*KK1NUS-1  ) 

I")  28  1 = 1 « !<TOP 

no  2-3  0 = 1.  NT  OP 

2h  U=U-  <4*i-2*K(nMlS-l>*  < 4*0-2 *NMll\lUS-l  ) * 

C PH  INN  ( 2*1-?  *k  hi  NUS-  1 . 2 *0-2  *NK  I IvlU.S-l  ) 

30  T=  I +U*iUh/2  ) * A ( m/2  ) 

w = w* ( 1+GOU) *AOH*KCUE*T 
T = 0 

00  AU  K = ST  EP  « S » S TFP 
no  AO  '\i  = SlEP.S.sTEP 
Nv)  I NUS=  C hi-  l ) - ( ( M-1  ) /?  ) *2 
(VTOP  = ( N-l  ♦ NillNllS  > /2 
U=-PK 1NU ( K »N ) 

00  30  1 = 1 .NT  OP 

38  U=U-  ( A * l-2*NMM|S-l  ) *PH  1 NO  ( K . 2*  I -MM  NUS  -1  ) 

AO  T=  r + U*h  ( K/2  ) ♦ A ( ,\i/2  ) *K  * ( K*1  ) *N 
P=u-GNU*AOH*KCOF*T 
T = 0 

00  30  K=b«EP.S,sTEP 
00  30  N = S IEP.S  »STEP 

U=-PMTAiU(K-l  «M*  1 )+PK  INN  (K-l  • N- 1 ) ♦PKINN  C K + l .N*l  ) -PKINN  ( 

K+l.N-1) 


MAUL'K  PWoOKArt  --  EVEN  (JUCALlNb  OF  A SPHEHILAL  SHtLL. 


fiPAKTY  = i\l-(M/2)*? 

NT0P=<N4MPAKTY>/2 
DO  H-)  1 = ] .NTOP 

4t>  U=U  + K*  <*\+l  )*(4«i-?*N|>a«TY-1  )*  ( -PK  INN  ( K- 1 «2»  I -NP  ARTY-1) 
C ■MJ«li\ll\i<  K4-1  ,<i«i-NPAK  f Y-l  n/(?»K*l  ) 

bU  T=T-U*AtK/2>*Al  m/2 I *K  *N* ( K + 1 >*(N*1)/<(?*K*1.0)*<  2*N+1. 

0 ) ) 

W = (rl*KHl)S*AOH*RCOE*  I 
T = 0 

DO  bl)  fsSlEP.S.STFP 

6U  T=r+P.*lM*(iM+l.  )*A(N/2)*H(N/2)/(2.*N+l.  ) 

W = W-M«COE*AOH*T 

W2=2-W 

Im3  = w2//> 

WHITE  Ib<  /n  ) U«M2»W3 

7U  F OK  MAT  ( 1 X . • D T ST  ANCE  B E T I*  F F M POLES  r »,Fe.5. 

C / • 1 X • ’TOTAL  POLAK  UEFLECTIQN  = »,F«.5t 

C /.lA.’SEMT  POLAH  REFLECTION  = »,F8.5> 

wKirf  (b,«(n 

fto  F OH  M A r ( lb ( bdX ,****•,/)  ) 

KE I UNN 
ENU 


■WW"  ' 


S a T C L L I It.  PttuttKAK  --  lNDbHALS  FOK  LLL  IKSOlU  PHObLEH  • 


C f'AIhi  HOOIlNt 

I • 1 i'IL  iv.  S lnlJ  Cl(4n.4  0),C2(40»40).C4(4ll>.S:6(40»4n>. 

C S / (HU  ) (S»Min  ,41?.  4 0 ) .S*»  (40 .40  ) »T1  ( API ) .72*  321  ) • 

C r 4(  VI ) .Th(3?1  >.Sb(4n.40).SlNE(4U> 
ht.K  = .e> 

hFEh>N  = Hfc.t  *Pt.f 

bNU=.3 

f.=HU 

1-11  = 4*/*  i 
MSMs2iM*hl- 

c sitp  ia:  coiMurL  imiE-Iihals  of  sin  u*xn. 
no  i in, him 

1 '11(1)50.0 

bU  2 llsl,  HI-1,? 

I =2uQ*I 1 
J5200-U 
TUI  >=2.n/Ti 

2 IK  J)  = -7  l ( I ) 

C SHIP  m:  COi.OlirF  iMTFbKALS  ok  kehnfl  * SIN  I i*xn 
C FC/P  A WIDE  KANbL  OK  VALUES  OF  1,  llSINb  SIKPSON'S  PULE. 
C IhE  VArtlA'SLI  JJJ  HKLOW  IS  4 WHEN  J IS  OOU, 

C 2 WHEN  J IS  M/ lM. 

Pl  = 4. 14  lb1*/ 

OA=Pl/lUOO 

00  b 1=1, HIM 

1 ?.  i 1 > = U 
T4( I )=0 

5 T 4 ( I ) = •) 

00  6 1 51,  ISM,? 

J5(-1)**( ( 1*4) /? ) 

T2(20u+I  )=j 
13(20041 )5j 
b T4(200*I)=J 
00  lO  J = 1 « 4**^ 

JJJ=242*( J-2*( J/2) ) 

Xi5j*0X 
Siv=Sl'v(  XI  ) 

CS=COS<  <1 ) 

0EI\l0*=l«EESO*CS4*24Sr'*4? 
bEN0H2  = UEN0i1*42 
UENUMsUEi\IOfl**l  ,b 

MEMO*'  45UKi\iO/4*l)rNOH 

00  10  I 15211. M*M, 2 
1=11-200 
SlsSIWI 1*A| > 

12(11)512(11  )4J.|J*SI/0tM0r3 
T4( 1 1 ) = I4( I I ) 4jjJ*SI /UEN0H4 
10  14(I1)=T4(I1) 4 J j«J*S  I /DENCF  2 
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y |\f|j|ij*n!5  • 


SATE  LI.  I It  I'KU'iH/i'l  --  INTI  UKALS  FOH  ELLIPSOID  PKOMLEM. 


ijo  12  1 1 = m 1 1 m * ? 

I 2(  ( l > = 1 ( I 1 > *0x*  . e>hfef,fet>7 

T 5 < [ 1 ) - M ( I I ) ♦i,y*.b‘>f,(,bb7 
T4(lI)=I4tIl)  *IJX*  . febfefefefe7 
1 = 11  - 2 0 U 
1 2 CONUNUt 

no  i<j  ii  = i. on 
1 =20  u - 1 1 

j=i?now  i 
12(I)=-TJ»(J) 

Til  1 ) =- I i ( J ) 

Id  I 4 ( I ) = - T 4 I J I 

C S1LP  2A:  CONPMTE  I pi  INTERNALS  SSII.J)  OIHLCTLY  • 

C 1 * S 1 N(,  SIMPSON »S  r>U|E.  Tut  VAKlAPLE  KKK  bELO* 

C IS  *4  WPEN  i\  IS  000,  2 LHEN  K IS  EVEN. 

110  40  1=1  .>1 
UO  '4  0 JS1.M 

40  ssi i ,j)=o.n 
00  AO  H=1 .494 
X I = k *i;x 

KXKs^4?*(X-|K/?1*?I 
SN=SIN(X1  ) 

CS  = UCS( Al  ) 

CSSN=CS*lS 

U E M 0 iv*  = H L E S Q * C S S n + S i4  ♦ S 14 
00  fell  J = ] .tf 
fell  SirjE  (J)=S1M(J*XT) 
r<ATIO=U.O 

00  HO  1 = 1,1" 

NAT1C  = KATI0*CS+C0S{  I 1-1  J*XI > 

00  ao  j = l 

«l)  S3  I I .0  )=ss  I I . J ) +KK»V*CSS0*HA1  10* SINE  ( J ) / Ut  NON 
00  4 U I=J  .1" 

00  4<)  0=  I • N » 2 

Ssl J « 1 >=SSI  I .0 ) *UX  * . fefefefeb67 
SSI  I « 0 ) =SS I J , I ) 

90  CONTINUE 

41  FUHhAT  (AX.iil  I?,i>  ) .E14.7) 

C STEP  2UJ  COPHurt  ImL  INTERNALS  S1-S4,  Sfe-Sl?. 

C A NO  T OE I •<  COMOlNA  I IONS  C1.C2.  A NO  C4,  AS  CESCIUUEO 
C IN  THE  SECTION  ON  EVALUATION  OF  INTEGRALS. 

DO  110  Jlsl.M 

00  111)  J=1,M 

1 = I I 420U 

Sl=-  I r 11  I4.J4l)-Tl  I 14.1-1  ) -1  1 ( I-J41  ) 41  1 ( I -J-l  ) ) /4 
S2=<T1  ( i+J4S)-3*Tl  ( 1 4J+1)  4,1*111  I4J-1  )-Tl|  1+J-i) 

C - T 1 I I-J  + s)4i-»Tl  | T-j-fl)-4*Tl(  I-J-l)4Tl  I I-J-3)  )/8 

S4=C  T1 ( I+J  + l )♦ ri  I I ♦ J- 1 ) 4T1 I I-J41 ) +T1 I I -J-l ) >/4 
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% 


(4 

'j 


Sudu  111  l,KUM<fl(1  --  INII  OMlS  FOK  ILL  ll’SOilJ  I'hOhLf.M, 


S4=-(  m 1 *'>+.}+  1 ) ♦!  1 ( T424J-1  >4  T1  ( ) + Tl  ( 1 42-J-l  ) 

C •i'«(Tl(l*JM)  + lt<[4J*)l  + Tlll*JU)+TlH-d*1l) 

C miI*?*J*l|MllI-?+J-l)  + THI-2-J+ll+U(l-2-lJ-l)l/ 

1 #> 

Ci(lI,J»=-J**24<bl4(  1-hLK  bQ  ) /HEF  f=U*S?  ) 

C 4J* ( S343* ( I-hLKSO ) /HFES(J*S4 1 

E -(HEFSG*Sb( I1»J)  ♦ GMU*Sl) 

110  C2Ul.ol:C1(ir.J)  ♦ 2*liMU*Sl 
no  120  ii=i.n 

i = I 1 ♦*011 

C ♦T'4(  I-2*l)+T4(  1-2-1  ) )/a 

SI  1 = - I T4  ( 1424*  1+3*  I1M  J*?*\  > 4T4  ( 1*2-1  ) ) 4T4  ( 1 + 2-3) 

c •?«ut(i  + M+i*(rmi*i)*T‘HHinT‘)(i-ii) 

C ♦ r«4  < T-24  3 )+3*  < T4( 1-2*1 I + T4 ( 1-2-1 ))  4T4 ( 1-2-3)  1/32 

S12=(T4(|*4*1)*t4(144-i)-4*(  r*4(  I + 2*  1 1 *T4  ( 1 ♦ 2 — 1 > ) 

C 4<i  * ( T4  ( I ♦!  )♦!“  <1-1 ) )-4*  ( r<?  U -241  >*T4  ( 1-2-1)  ) 

C 4T4( 1-4+1 )*T4< I-4-l ) )/3? 

120  C3( 11  )s-.b*dfcrSn*Sl0+BFLSt*(3-.5*GNU)*Sll4( l-.b*GNU42* 

liLESU  1 *S12 

no  i Ao  11  = 1 * ii 
1=11 *200 
no  130  jsi.m 
no  12b  * = 1 .11 

1 2b  SH ( I l . J.K )=- ( T2< I 4 J ♦ K 4 1 14 » ? ( ) IJ1K-1I -T?  < 14J-K41  )-T2(  14 

J-K-l  ) 

C -T2<1-J4K41)-T2( I-J4K-1 )4T2( t - J-K 4 1 1 41 p ( J - j-K- 1 

) ) /« 

Sbll I . J»=-( T2( 14J41 i-xp( i+j-1 )-12< I-J41 )4T2( I-J-1 1 >/4 
S'*  < 1 l . J>  = < Ti(  1 4J43  )-3*T3<  14J41  )43*T3<  I4J  — 1 1 -T3(  14J-3) 

C -T3<  r-,j43>43*T3(I-J4i)-3*r3(l-J-l)473(I-J-3)  ) 

/IF. 

130  S7(lI)  = -<I.Ml424l>4i.M|4?-l  ) -2*  (T3(14l)4T3(l-l)l 
C.  4 rb(  1-241  )*T3<  i-2-l  > )/« 

1)0  200  1 = 1.11 

200  UK  1 I L (0,201)  I ,C3< I ) ,S7< I » 

201  FOXMAT  ( 1 2 . t 14 . 7 ,E 14 . 7 . *>0X ) 

00  2 in  l = l,.l 

t<0  210  Jrl.rt 

210  4l<llt  (0,202)  l , J.C1  ( I , J ) , (2  ( I , J ) , SO  ( I , J ) » S9  ( I • J ) 

202  F OK ft AT  ( 12,12,t14,7,Ll4.7,E14.7.F14.7,?0X) 

kHITF.  ( H ,203  ) ( ( ( ( 1 ,J,K«SM(I«J,K)  1,1  = 1, J),J  = 1,K)  , K = 1 • M 

) 

2U3  F0KKAT(H(12,12.T2,E14.7I  ) 

S f IJK 

end 


Tn 


mi 
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SuTLL.l  l IK 


1N1L0HALS  F OK  SPHM<t  PKOHLLM 


C MAIM  Pool  INL 

iml<’i:i<  nirsi , p h i<  l * i * a r<  r i . s t i p,ss,ssh.ssm 
liOUuL  t PptClSION  U.v/,rt(?/*0) 

LO  Jt'STzi 
H lti>  T = S 1 
STKP  = L 
WA=2*H ILS I 
iW  April  A + 1 
SS  = 4 *H  Its  T 
SSP  = SS+  1 
SSM=SS-1 
A ( M A ) = 1 

Du  go  i=nap.ss 

J=1  -C.A 
A ( J ) =>l 

fell  A ( 1 ) = ( P*  I-SSP  ) * A ( I -1  ) / ( 1 -HA  ) 

A ( HA ) =1 
h APPAzO 

DO  PliO  LzLOvJLS  f ,HILST  tSTF  P 
DO  POO  rizGl  K P « I .Site 
DO  POO  HzSTLP  , t-l , ST  Lf 
P/W<L  = L-  < l./P  ) 

LA=(L4pAKL)/P 

il  < (L  + ii+i\.)  .c-j.  ( (i  ♦rHN>/2)*;n  go  tu  poo 

IK  I L . G r , ( 1 + I-.  ) ) GO  10  POO 
P A M (J  = IM  - ( K / ? ) -»  p 
MAz ( w - P A K O 1 /P 
V = U 

DO  IMi  1=1  t LA 
1 A z A 4 [-1-P*PARL 
1 B z P * I -1-PAPL 
I Czil-1 

ill=(  IO*IC*M)/?  + mA 
I L = H ♦ 1 

1F=< iH+ir+N) /?4MA 

u=o 

IK  (liA.LO.O)  GO  TO  150 
DO  14IJ  J=  1 ( DA 

ja=a  * j-5+p  *paph 

Ut'  = P*  J-P  + PAKi  J 
JC=( lP  + iOJH) /P+HA 
JUz  ( Ih+  it  + Jti  ) / 2 + MA 

140  U=U-»-JA*P*  ( A T JO-Ib  ) *A  ( JC-iC  ) *A  ( JC-Jb  » / ( ( 2*JC-SSM  )*A  ( jC  ) 

) - 

C A ( Jll-I'-n  ♦ A ( jD-  if  ) * A C JO-JQ  ) / ( ( ?.*  JD-SSM)  ♦ A ( j|t  > ) ) 
lbO  V = V+1A*  A ( JO-IH)  *A  ( TD-IC  )*A  ( 1U-IV)/  ( < 2*  IU-SSM  ) * A ( I 

0 ) ) - 

C AUK-lB)*A(IK-It)»A<IF-N)/<  ( 2*  I F -SSf'; ) * A ( IF)  ) ) ♦U  > 


I 


SMlLLlIlt  HR()G*A»1  --  IMlGKALS  FOR  SHHtKt  FKObLLM. 


V= )*W/  > 

K ApHA=rt  AfHA-fl 
UK  lit  (o.H'*)  L.n.N.V 
V5S  FjlblAT  !UilMP.hl4,7) 

200  CONTlNbt 

WKITK  lfc,201>  KflF'PA 

201  FOKi-./U  (POX. IF,, • RLCOKOS  WKHIEImM 
iU  P 

cNJ 


COPY  mu$f T"  " 


nr* 
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SAltLLiffc  HKUGKiVM  --  H.KlulAl\i  FOK  SPntKt  PBuHLt-M. 


C M A 1 U KOUTiMt 

Ulr'lLtN.biOIJ  XHuSlinO)  , YPOS(IOO)  ,A(bO>  .Hb(bU) 

IwltGiP  S.SfF.P 
KtAL  NL 

UP  ( IM  , x ) =.im*  (P  i M-  1 . X ) -x*P(\  , * ) > / ( ( 1-X*X  ) **.b  ) 

UDP  ( N«  X > =0*  < < -1  -IM+l4*V  *x  ) *P  < IV  t X ) 4X»P<N-1  i X ) > / < X*X-1  ) 
i<FAU  ( b < 10  ) AOH.MJ.PHOS.X  . XX 

s = > 

s tkp=xa 

H t".  A U (b.lO)  ( A ( T ) . I = STl'PtS.SUP) 
ktAU  (b.lUl  (PS( 1) . 1=STLP«S,STLP) 

III  PUB  I'i  A T ( t 17. 10  .a  AX  > 

KCOL-  ('‘./l  A Oil*  *?*f«UUT  ) ) * ( 1 + AU/  ( 2*AUH*KOCT  ) ) 

F = KC0F  * A ( i H 
Pl=3.1AlbyA 
L = 4 

C)X  = PI  / ( 1 (i II  *L  ) 

LL  = ll)U«L-l 

*'0  A — 0 

UO  ion  1=1 iLL 

x 1=1 *nx 
Y=COS  < XI  ) 

tit  r a= o 

F’SIS  = U 
PSISP=U 

L'O  yu  IS^SIPP.S.STCP 
U.NP=DP <N* Y ) 

Btl  A = mP  T A + A ( f.< ) * TtMP 
PSIS  = PSlS  + f«b  (f.  ) *TE.MP 
SMJ  PS  I SP=PSI  sP*.JS  ( M ) *UUP  ( l\| « Y ) 

X = S1A  ( A 1 > 

Y =-  Y 
SX  l = X 
CX  I = -Y 

T£MP  = -Bt.l  A*CXT-{  1-MU  ) /2  *KHOS  *F*SX  I + F *PS  1 S*t  X I -MJ*F  *PS  I 

SP*SX 1 

C ♦ ( 1-fvil) ) / ?*KPO«?*F*PK'l  A*CX  T-r-IU*HHOS*P*RtTA*SXI**?*CXl 
C -,b*OtTA**?*axI-F  *l>SlS*HLTA*CXI**?/SXl*MU*F*PSISP*bt 

r A*CXI 

C ♦ WU*HHni>i‘F*Pfc.TA**2*SX  I*CXl**2*F*PSIS*fltTA*SXI 
TFKP=OX*7tMP 
MJAswUA+ILN.H 

IF  ( I r.  ( I/L  )*l  ) 00  TO  1 00 

XPOS<  l/L  I =X*F*  ( (NU-1  ) /2*HHOS*SX I ♦PS1SP*S>  1 - rJU*PR  I S *C  X I 

♦ I A 

C *SX  T ♦ ♦2’*Cx  I ) 

YP06<  I/L  >=Y*W()A-1  FOP/2 
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END 

DATE 

FILMED 


NL 


77 


SmULLIU  PKOGHA*  OLKlUIAb  FOH  SPHfcKt  KKOULtrt. 


loo  titan  hut 

1 POS  ( 1 OU  I si  ♦*iOA 
AO  JUJ»T=  ( 1 -VPOR 1 1 00 ))  /2 
oo  lub  lri. inn 
105  YPOSm  = tPoSU  UaOJUST 
DEF  Li'C  s * * 1 1 - 1 POS  1 1 0 0 ) ) 
olST=?*rHO«?uno» 

WHllt  (o.lOOO)  nlSl.l'tFUC 

1000  FOHrtAT  1 1 A • •ulST  »• .€i4.7./.lXi*0EFLtC  *«.E14,7» 
STOP 
ENO 


Sm11.LI.uk  HHUUHAf* 


KUO*  srtiLKC  TO  tLLlHSOlO. 


■ 


I 


C hAlU  HoUTiNt 

C A SUUTIUN  Fn«  04  T A*  EXPRESSED  IKS  TERMS  OF 
C LEGtAOMt  MINCT I ntoS « IS  TRANSLATED  INTO 

C A SOLUTION  FOR  OtTA,  EXPRESSED  IN  TERMS  OF  SlNtS 

INTEGER  SStSTEP.SNEto 
Ml AL  MitLAOH 
DIMENSION  ? ( ?0  ) ( A ( 20  ) 

|)P|N,k  )=N*(P(N-1  tX)-y*P(N,x> )/<  <1-X*»  )»*.5) 

UO  1 1=1, ?0 

1 ^uuo.lt 

HF.AU  (St?)  NHOST tNUt AOMtSStSlEPtSNLW 

2 FORMAT  (4EIS.7.M?) 

READ  (St3)  (A(l)tlsSTEPtSStSTEP) 

^ FORMA!  (fclS.7) 

sor=is«(i-No»*? » )**.s 

LAOrt=)4WU/(  (2*A')H)*S0T  ) 
RHoCRf=2*LA0H/(s0T*A0H*A0H) 

KHOSHNUSTAHHOCHI 

x=swtb 


HLE  = l 

A ITCHsl'LF./ AOM 

WH  lit  l*t2M)  V«(>HOtNUtHtF  t AITCH 
DX=3.)4lbV3/ir>0 
10  UO  -?«  M=ltS*) 

A=H*IJX 

XlsCOS(A) 

UU  2l»  N=S IFPtSS.STEP 
Y=OP(NtX| ) 

DO  20  1 = 1 tSWl'W 
SINE=SIN(1*X) 

20  2U  )=?(  i (♦A(N)*Y*SlNt 
DO  50  1=1 t SMEW 
A(l)=?(i)«2*0X/^.lHlb0i 
WHITE  (*•  = *»>  Alf) 

20  FORMA?  I FI**.  7) 

30  CONTINUE 
STOP 
ENL 


COPY 

PERMIT  Fii 


iBLE  TO  000  MB  K8 
V LEuEil  i mXtiA 
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S/>  TtLLl  It  PKON»w*rt 


FHOM  SPHERE  TU  tLLlrSOlb 


K*NC I 10U  PCN.X ) 

C TM*  N-IH  Lt-bcM'jKF  pOltUOMIAL  iVAlUAICu  A1  X. 
Hal 

IP  (K.LO.l)  P=X 
i»-  (fj.Lr.«;>  hetiihw 
s=i 
P=X 

DO  1U  I=?.N 

H = S 

S=l> 

10  P=ll2tI-ll*X»S-(I*ll*R)/l 

KETUHA, 

END 


i 

I 

I 
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